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Abstract 

Initial  efforts  to  characterize  the  scattering  dynamics  of  B  +  H2  focus  on 
computing  scattering  matrix  elements  for  the  fine  structure  transition  B  (  P1/2 ) 

B  ( 2P3/2  )  in  collisions  with  H2,  allowing  for  rotational  excitation.  Using  a  new 
application  of  the  time  dependent  Channel  Packet  Method  (CPM),  reactant  and  product 
wave  packets  are  prepared  in  the  asymptotic  limit  on  the  B  (  P1/2 )  and  B  (  P3/2 ) 
surfaces.  They  are  propagated  using  the  split  operator  method  together  with  a  unitary 
transformation  between  the  diabatic  and  adiabatic  representations.  Scattering  matrix 
elements  are  computed  from  the  Fourier  transform  of  the  correlation  function  between  the 
evolving  wave  packets.  These  computations  directly  support  the  Air  Force  Office  of 
Scientific  Research  (AFOSR)  Molecular  Dynamics  program  and  the  Air  Force  Research 
Laboratory  (AFRL)  High  Energy  Density  Matter  (HEDM)  program.  In  particular,  the 
CPM  is  well  suited  to  handle  non-adiabatic  molecular  reaction  dynamics  on  multiple 
potential  energy  surfaces,  as  encountered  in  the  dynamics  of  a  wide  variety  of  molecular 
systems,  including  B  +  H2 .  Further  motivation  for  investigating  the  specific  dynamics  of 
B  +  H2  stems  from  the  potential  application  of  solid  molecular  hydrogen,  doped  with 
boron  atoms,  as  a  high  energy  rocket  propellant. 
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INELASTIC  SCATTERING  MATRIX 


ELEMENTS  FOR  THE  COLLISION 
B(2P1/2)  +  H2(y)  ->  B  ( 2P3/2 )  +  H2  (f  ) 


I.  Introduction 


Motivation 


The  Air  Force  High  Energy  Density  Matter  (HEDM)  program  investigates  the 
development  of  improved  propulsion  and  explosive  systems.  A  principal  goal  in 
propulsion  is  to  improve  the  specific  impulse  ( Isp )  of  the  propellant.  Efforts  to 
synthesize  a  novel,  high  performance  rocket  propellant  using  boron  atoms  trapped  in  a 
matrix  of  solid  molecular  hydrogen  [1]  have  recently  generated  a  great  deal  of  interest. 
Theoretical  calculations  have  shown  that  if  one  could  achieve  doping  of  5%  molar  atomic 
concentration  of  boron,  the  specific  impulse  would  show  a  21%  improvement  over 
standard  liquid  oxygen  -  liquid  hydrogen  propellants  [1].  Thus  there  is  great  interest  in 
the  trapping  of  boron  and  other  atomic  species  in  solid  molecular  hydrogen  matrices.  As 
of  1994,  however,  only  weak  B  atom  molar  concentrations  (0.01%)  could  be 
achieved  [1],  Hence  the  spectroscopy  and  reactive  dynamics  of  systems  of  boron  and 
hydrogen  are  of  great  interest.  The  particular  scattering  event  we  study,  motivated  by  this 
interest,  is  the  collision  B  (  2Pi/2  )  +  H2  (j )  -»  B  ( 2P3/2 )  +  H2  {f  )  .  This  collision  is 
not  reactive,  but  it  is  inelastic  due  to  the  possibility  of  the  fine  structure  transition 
B  (  2Pi/2  )  — »  B  ( 2P3/2 )  and  the  possibility  of  the  rotational  transition  H2  (  j  )  -> 

H i(f). 
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The  dynamics  of  a  collision  between  a  single  boron  atom  B  and  a  hydrogen 
molecule  H2  are  completely  described  by  matrix  elements  of  the  scattering  operator 
(S-matrix  elements).  The  absolute  values  squared  of  S-matrix  elements  yield  the 
probability  of  transition  or  reaction  from  a  given  initial  state  to  a  given  final  state. 
Together  with  the  distribution  of  reactant  states,  S-matrix  elements  can  be  used  to 
compute  cross  sections  and  reaction  rate  constants.  These  quantities  are  expected  to 
contribute  to  large-scale  synthesis  efforts  of  boron  in  hydrogen  matrices. 

Several  authors  have  treated  analogous  systems  of  atoms  making  fine  structure 
transitions  during  collisions  with  diatomic  molecules  or  noble  gases.  Chu  and 
Dalgamo  [2]  studied  fine  structure  transitions  in  collisions  of  C+  (2P)  and  H2,  where  the 
complex  was  assumed  to  be  collinear  and  the  hydrogen  molecule  was  constrained  to  its 
ground  rotational  state.  Flower  and  Launay  [3,4]  extended  the  study  of  C+  +  H2,  lifting 
the  collinearity  restriction,  and  considering  the  possibility  of  simultaneous  rotational 
excitation  of  the  H2  molecule.  Also  related  are  several  studies  of  fine  structure  transitions 
between  alkali  atoms  and  inert-gas  atoms  including  Na  (  Pj)  +  He  [5],  and  Na  (  Pj)  + 

Ar  [6].  Rebentrost  and  Lester  [7-10]  carried  out  a  study  of  the  F  (2P)  +  H2  system, 
including  the  construction  of  the  potential  energy  surfaces,  transformation  to  an 
approximate  diabatic  basis,  and  a  time  independent  calculation  of  scattering  matrix 
elements  and  cross  sections  for  the  fine  structure  transition,  including  rotational 
excitation  of  the  hydrogen  molecule.  In  this  system,  the  spin-orbit  splitting  of  fluorine 
( 404  cm'1 )  [9,10]  roughly  corresponds  to  the  energy  required  to  make  ay  =  0  -»y  =  2 
rotational  transition,  so  there  was  a  near  resonant  effect  that  enhanced  that  cross  section. 
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Problem  Statement 


To  make  a  propellant  that  achieves  the  theoretical  improvement  in  specific 
impulse  described  above,  one  must  be  able  to  achieve  higher  concentrations  of  B  atoms 
trapped  in  the  H2  matrix.  To  do  this,  an  understanding  of  the  dynamics  of  B  atoms  in  the 
H2  matrix  is  important.  A  first  step  in  studying  these  dynamics  is  to  study  the  dynamics 
of  a  single  B  +  H2  collision,  characterized  by  scattering  matrix  (S-matrix)  elements.  We 
compute  S-matrix  elements  as  a  function  of  energy  for  the  fine  structure  transition 
B  (  2Pi/2  )  ->  B  ( 2P3/2 )  during  collisions  with  H2,  while  allowing  for  the  possibility  of 
rotational,  but  not  vibrational,  excitation  of  the  hydrogen  molecule. 


Approach 

S-matrix  elements  can  be  computed  by  a  variety  of  methods.  We  will  use  the 
time-dependent  Channel  Packet  Method  (CPM).  Time  dependent  methods  give  S-matrix 
elements  as  a  function  of  energy,  for  one  set  of  quantum  numbers  that  identify  specific 
initial  and  final  states.  The  time  dependent  Channel  Packet  Method  uses  propagation  of 
wave  packets  that  explore  the  interaction  potential  to  derive  S-matrix  elements  as  a 
function  of  energy  for  the  specific  initial  and  final  states.  Previous  calculations  for  similar 
reactions  have  used  time  independent  methods  such  as  close-coupled  methods  [11-13]. 
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To  compute  S-matrix  elements  one  must  know  the  interaction  potential. 
The  potential  energy  surfaces  governing  the  dynamics  of  the  B  +  H2  interaction  have  . 
been  computed  by  Alexander  [14].  They  display  a  strong  barrier  to  formation  [14]  of  the 
BH2  molecule,  and  at  low  energies,  the  formation  of  a  van  der  Waals  complex,  B  -  ■  -Kb  ,  is 
possible  [15].  The  interaction  of  B  and  H2  is  complicated  by  the  fact  that  the  ground  state 
of  B  ( 2s22p  2P  )  is  sixfold  degenerate,  if  spin-orbit  coupling  is  neglected.  This 
degeneracy  means  that  the  interaction  takes  place  on  multiple  coupled  potential  energy 
surfaces,  rather  than  a  single  potential  energy  surface,  adding  to  the  computational 
difficulty.  The  interaction  dynamics  are  described  by  three  potential  energy  surfaces 
(independent  of  spin)  corresponding  to  the  three  2p  states  of  the  boron  atom  as  it  interacts 
with  the  hydrogen  molecule.  The  six  total  surfaces  describing  the  interaction  of  B  (  Pj) 
withH2  are  only  separated  by  a  constant  spin-orbit  splitting  term  in  the  limit  of  infinite 
separation  of  the  boron  and  hydrogen,  but  they  are  split  at  shorter  range  by  the 
interaction. 

Section  II  gives  a  brief  background  on  the  theory  of  S-matrix  elements,  the  use  of 
the  channel  packet  method,  and  the  concept  of  adiabatic  and  diabatic  potential  energy 
surfaces  to  model  the  interaction.  In  Section  III,  we  identify  the  coordinate  systems  and 
the  scattering  Hamiltonian.  The  Hamiltonian  is  then  represented  in  a  complete  basis.  We 
then  make  simplifying  approximations  by  neglecting  the  vibrational  degree  of  freedom  of 
the  hydrogen  molecule,  limiting  the  highest  rotational  state  accessible  to  the  hydrogen, 
and  using  the  “Centrifugal  Sudden”  approximation.  Section  IV  describes  the 
implementation  of  the  channel  packet  method,  and  discusses  specific  numerical  issues, 
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including  wave  packets,  the  numerical  fidelity  of  the  calculation,  and  methods  to  validate 
the  calculated  S-matrix  elements.  Section  V  presents  the  results  of  the  scattering 
calculations.  Results  are  first  discussed  for  the  simplest  case:  where  the  hydrogen  is 
restricted  to  its  ground  (j  =  0  )  rotational  state.  S-matrix  elements  where  rotational 
excitation  is  allowed  are  then  presented. 
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II.  Background 


S-Matrix  Elements 


A  scattering  problem  can  be  characterized  by  the  scattering  operator  S ,  which 
gives  the  asymptotic  ( t  -»  +°o)  product  state  |  xVproduc^  in  terms  of  the  asymptotic 

( t  -oo  )  reactant  state  \^reaclan)  by 


\JJ  \  _.  c  vp 

A  product  j  u  L  reactant 


(1) 


Knowledge  of  S  gives  one  all  information  of  experimental  interest  [16],  since  many 
reactions  typically  take  place  on  very  short  time  scales  and  the  intermediate  states  are  not 
experimentally  observable.  Matrix  elements  of  the  scattering  operator  can  be  used  to 
obtain  the  probability  of  a  final  state  given  the  initial  state: 


(2) 


We  desire  scattering  matrix  (S-matrix)  elements  where  the  initial  and  final  states  are 
asymptotic  momentum  eigenstates,  labeled  by  the  reactant  and  product  momenta  &rand 
channel  quantum  numbers  y 


6 


(3) 
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r'r 

±kr<±kr 


(k'r  1  s  |  k  r) 


2 

The  amplitude  squared  |  S  |  of  the  S-matrix  elements  gives  the  probability  that  the 
transition  will  take  place.  From  S-matrix  elements,  one  can  compute  scattering  cross 
sections  and  the  corresponding  rate  constants. 


The  Channel  Packet  Method 

The  time  dependent  channel  packet  method  (CPM)  [17,18]  is  based  on  the  Moller 
operator  formulation  of  scattering  theory  [16].  The  Moller  operators  are  defined  by: 

=  lim  exp  (^-)  exp  (  )  ( 4 ) 

r->+«  fi  h 

where  H  is  the  Hamiltonian,  and  Hr0  is  the  asymptotic  channel  Hamiltonian  labeled  by 
the  quantum  numbers  contained  in  y  .  The  Moller  operators  have  the  property  that 
S  =  Q_  t  Q+  .  The  Moller  operators  can  be  used  to  define  Moller  states  in  terms  of 
initial  reactant  and  product  states,  defined  arbitrarily  at  t  =  0: 


vF/\  =  Qr  x\>r  \ 

T±  /  —  1  in,  out  / 


(5) 
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The  initial  reactant  (product)  states 


xur  ' 

T  in, out  i 


are  typically  chosen  to  be  Gaussian  wave 


packets  entering  (leaving)  the  interaction  region,  so  that  they  possess  momentum 
components  that  will  give  S-matrix  elements  (Eq.  (  3  ) )  that  are  of  interest.  Once  the 
Moller  states  are  constructed,  the  correlation  function  Cyr(t)  is  then  computed  as  one  of 
the  Moller  states  is  evolved  under  the  full  Hamiltonian: 


C„(<)=(^r|exp(-^p)|'Fl) 


(6) 


The  Fourier  transform  of  the  correlation  function,  appropriately  normalized,  will  give  the 
S-matrix  elements  as  a  function  of  energy  [19]: 


<±ky 


{My  Myf2  V-(.±k' r.)  ?l+(±kr) 


exp  (~~~)  Cyy{t) 


(7) 


In  Eq.  ( 7  )  E  is  the  total  energy,  /labels  each  reaction  channel,  juy  and  //Y  are  reduced 
masses,  and  //_  and  r]+  are  momentum  space  expansion  coefficients  used  to  define  the 
initial  reactant  and  product  wave  packets.  The  initial  wave  packets  must  be  prepared  to 
have  either  all  negative  or  all  positive  momentum  components  for  the  above  formula  to 
be  applicable.  A  more  general  formula  [17]  is  required  when  this  is  not  the  case. 
Typically  it  is  most  convenient  to  prepare  the  initial  wave  packets  as  described  above. 
However,  one  difficulty  with  the  channel  packet  method  occurs  for  low  energy,  where  the 
initial  packets  do  not  possess  significantly  large  momentum  components.  At  these  low 
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energies,  the  numerical  result  becomes  invalid,  because  one  is  dividing  a  small  number 
(the  Fourier  transform  of  the  correlation  function)  by  another  small  number  (the 
expansion  coefficients).  The  use  of  wave  packets  with  both  positive  and  negative 
momentum  components  may  help  with  this  situation,  but  it  complicates  the  evaluation  of 
the  S-matrix. 

To  apply  the  channel  packet  method  to  a  scattering  problem,  one  must  define 
initial  reactant  and  product  wave  packets  that  contain  momentum  components  that 
correspond  to  energies  for  which  the  S-matrix  elements  are  desired.  In  our  multiple 
surface  problem,  the  initial  reactant  packet  can  be  defined  on  one  surface  and  the  initial 
product  packet  on  the  other  surface,  so  that  the  transition  probability  between  the 
surfaces,  which  represent  different  initial  and  final  sets  of  quantum  numbers,  can  be 
found.  Next,  the  reactant  and  product  Moller  states  are  numerically  computed  by  the 
action  of  the  Moller  operators.  The  Moller  operators  first  propagate  the  reactants 
(products)  backward  (forward)  in  time,  away  from  the  interaction  region,  under  the 

asymptotic  channel  Hamiltonian  H£  .  This  can  be  done  analytically  if  the  initial  states 

are  chosen  to  be  Gaussian  wave  packets,  which  is  often  the  case,  for  simplicity.  The 
analytic  propagation  is  to  a  time  x  that  is  large  enough  so  that  the  wave  packet  is  then 

(numerically)  in  the  asymptotic  limit,  where  the  full  Hamiltonian  H  is  equivalent  to  the 

asymptotic  Hamiltonian  HI  .  Then,  the  wave  packet  is  numerically  propagated  in  the 

* 

opposite  direction  of  time  under  the  action  of  the  full  Hamiltonian  H ,  back  to  t  =  0, 
where  it  is  then  the  Moller  state. 
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Once  the  Moller  states  are  created,  the  next  step  is  to  compute  the  correlation 
function  between  them  while  one  of  them,  say  the  reactant  Moller  state,  is  evolved 
backward  and  forward  in  time  under  the  action  of  the  full  Hamiltonian.  This  is  done  with 
the  same  numerical  propagation  method  used  for  the  computation  of  the  Moller  states.  A 
time  range  [  zi ,  r+  ]  is  chosen  to  approximate  the  true  range  (-00,  co)  by  stopping  the 
evaluation  of  the  correlation  function  when  it  decays  to  numerical  insignificance,  when 
the  evolving  reactant  Moller  state  has  finally  completely  exited  the  interaction  region. 

After  the  full  correlation  function  is  computed,  it  can  then  be  used  in  Eq.  (  7  )  to 
compute  the  S-matrix  elements.  To  evaluate  the  S-matrix  elements  as  a  function  of  the 
total  energy  of  the  system  we  require  the  implicit  dependence  of  the  expansion 
coefficients  rj  (k)  on  energy,  given  by 


kr(E)  = 


h 


,1/2 


(E~Er) 


(8) 


where  Er  is  the  asymptotic  energy  of  the  channel  labeled  by  y,  for  example  the  rotational 
and  electronic  fine-structure  energies  in  our  problem. 

Previously,  the  CPM  has  been  used  in  scattering  problems  involving  rovibronic 
transitions  on  a  single  electronic  surface,  including  computation  of  state  to  state  S-matrix 
elements  for  the  collinear  H  +  H2(  v)  H  +  H2(  ft)  reaction  [19-21].  Another  two 
dimensional  calculation  computes  S-matrix  elements  for  OC  +  OH  (  v=  0) 
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OCO  ( v-  0)  +  H  [22].  In  a  three-dimensional  calculation,  Dai  and  Zhang  [23]  have 
used  the  CPM  to  compute  S-matrix  elements  for  the  H  +  O2  reaction.  The  CPM  has  also 
been  used  to  develop  approximate  methods  for  computing  S-matrix  elements,  including  a 
semiclassical  method  [24]  and  a  method  adapted  to  work  with  a  time-dependent  self 
consistent  field  approach  [20]. 


Atomic  Units 


Atomic  units  [25]  are  convenient  for  atomic  and  molecular  calculations.  One 
atomic  unit  of  mass  is  the  mass  of  the  electron,  one  atomic  unit  of  energy  is  twice  the 
ionization  energy  of  the  hydrogen  atom,  one  atomic  unit  of  distance  is  the  Bohr  radius  ao , 
and  one  atomic  unit  of  time  is  the  time  for  an  electron  in  the  first  classical  Bohr  orbit  to 
travel  one  atomic  unit  of  distance.  The  atomic  unit  of  angular  momentum  is  fi ,  which 
conveniently  allows  fi  to  be  omitted  from  equations  presented  in  atomic  units. 


Adiabatic  and  Diabatic  Potential  Surfaces 


To  solve  a  scattering  problem  one  needs  a  reliable  description  of  the  potential 
energy  surfaces  that  describe  the  interaction  between  the  scattering  partners.  A  basic 
approximation  required  to  even  discuss  the  concept  of  a  potential  energy  surface  is  the 
Bom-Oppenheimer  approximation  [26],  which  allows  one  to  use  potential  energy 
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surfaces  to  describe  the  nuclear  motion  under  certain  assumptions.  In  this  approximation, 
the  nuclear  position  operators  are  demoted  to  parameters  and  the  electronic  Schrodinger 
equation  is  solved  for  the  electronic  wavefunctions  and  energy  levels,  as  functions  of  the 
nuclear  positions.  The  electronic  energy  levels  as  a  function  of  nuclear  position  form 
adiabatic  potential  energy  surfaces  and  define  an  adiabatic  basis. 

Using  this  approximation  involves  neglecting  terms  that  take  the  form  of  matrix 
elements  of  the  nuclear  momentum  operators  in  the  full  Hamiltonian  between  the 
electronic  eigenstates.  This  means  that  the  nuclear  kinetic  energy  operator  in  the  full 
Hamiltonian  is  not  diagonal  in  the  adiabatic  representation.  The  approximation  is  good 
when  the  energy  levels  are  well  separated  compared  to  the  energy  of  the  nuclear  motion. 
This  can  be  a  good  approximation  for  interactions  with  closed-shell  atoms,  in  which  the 
electronic  energy  levels  are  well  separated,  on  the  order  of  eV.  However,  for  interactions 
with  open-shell  atoms  such  as  B  (2P) ,  the  2Pj  level  is  nearly  sixfold  degenerate,  since  the 
energies  for  the  2Pi/2  (twofold  degenerate)  and  2P3/2  (fourfold  degenerate)  levels  are  close 
(the  spin-orbit  splitting  constant  £,  for  B  is  10.7  cm'1)  [14,27].  Thus  the  off-diagonal 
matrix  elements  of  the  nuclear  kinetic  energy  operator  cannot  be  ignored.  They  can  be 
treated  approximately  by  transforming  to  a  diabatic  representation  [26,28-32],  where  the 
nuclear  kinetic  energy  is  diagonal  but  the  potential  energy  is  no  longer  diagonal.  The 
transformation  to  the  diabatic  representation  can  be  derived  by  solving  a  differential 
equation  involving  the  above  mentioned  matrix  elements  [28]. 
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The  determination  of  the  potential  energy  surfaces  for  molecular  interactions  is  a 
challenging  problem.  One  may  use  ab  initio  methods  to  solve  the  electronic  Schrodinger 
equation  at  various  nuclear  locations,  obtaining  the  electronic  eigenvalues  and 
wavefunctions.  Examples  of  these  calculations  are  self-consistent  field  methods  and 
configuration-interaction  methods  [11],  which  are  both  used  to  solve  the  Hartree-Fock 
equations.  However,  certain  regions  of  the  surfaces  may  be  known  experimentally  to 
better  accuracy  than  is  available  theoretically,  so  another  approach  is  to  use  experimental 
data  to  construct  analytic  models  for  which  the  potential  matches  the  experimental  values 
at  the  known  locations. 

To  construct  the  (approximate)  diabatic  representation,  one  must  know  the 
electronic  wavefunctions,  which  implies  that  one  has  already  calculated  the  adiabatic 
energy  levels.  Thus  in  an  ab  initio  treatment  of  the  potential  surfaces  for  a  problem  such 
as  ours,  it  is  essential  for  a  representation  of  the  diabatic  potentials  to  be  given  with  the 
adiabatic  potentials,  or  one  would  have  to  redo  the  entire  potential  surface  calculation  just 
to  obtain  the  transformation  from  the  adiabatic  to  diabatic  basis.  Without  the 
transformation,  there  is  not  enough  information  to  model  the  interaction  fully. 

The  Hamiltonian  matrix  in  the  diabatic  representation  will  be  obtained  in 
Section  III.  The  potential  energy,  including  the  angular  effective  potentials,  is  not 
diagonal  in  the  diabatic  representation.  The  diabatic  potential  matrix  Vn  (R)  is 
diagonalized  to  obtain  the  adiabatic  representation,  which  is  used  in  the  wave  packet 
propagation.  The  unitary  transformation  U  (R)  that  accomplishes  this  diagonalization  is 
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obtained  by  placing  the  eigenvectors  of  the  matrix  in  columns  of  a  new  matrix,  and  can 
be  used  to  transform  the  wavefunction  between  the  adiabatic  and  diabatic  representations: 


VD  VoiR)  =  UlAVoUlA  yD(R) 

=  U  VA  Ut  VoiR)  (9) 

=  U  VA  yA(R) 

where  ^d(R)  and  ^(7?)  are  the  diabatic  and  adiabatic  nuclear  wavefunctions  (column 
vectors)  and  VA  (R)  is  the  diagonal  adiabatic  potential  matrix.  Hence  the  potential 
energy  Vd  (or  any  power  of  VD ,  and  hence  any  function  of  Vd)  acting  on  the  nuclear 
wavefunction  in  the  diabatic  representation  T^i?)  can  be  evaluated  with  Eq.  ( 9 )  by 
applying  the  Hermitian  conjugate  of  the  unitary  transformation  U  (R)  to  the  state, 
applying  the  diagonal  adiabatic  potential  energy  VA  (R)  operator  and  then  transforming 
back  with  U  (R) .  The  benefit  of  this  transformation  approach  is  that  during  the  wave 
packet  propagation,  the  Hamiltonian  is  exponentiated,  which  is  straightforward  for  a 
diagonal  operator.  The  unitary  transformation  that  diagonalizes  the  diabatic  potential 
energy  matrix  will  be  a  function  of  the  nuclear  coordinates.  However,  any  unitary 
transformations  that  are  not  functions  of  the  nuclear  coordinates  may  be  applied  to  a 
diabatic  representation  and  it  will  remain  diabatic  (the  nuclear  kinetic  energy  will  still  be 
approximately  diagonal)  [26]. 
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III.  Representation  of  the  Hamiltonian 


Overview 


To  treat  the  scattering  problem,  we  must  identify  a  coordinate  system,  the 
Hamiltonian,  and  a  basis  with  which  the  Hamiltonian  is  represented.  Much  of  the  general 
procedure  for  atom-atom,  atom-diatom  and  atom-molecule  collisions  is  developed  in 
several  papers  and  books  [1 1,33-47].  We  outline  the  key  steps  here  as  applied  to  the 
B  +  H2  scattering  problem.  Most  approaches  are  designed  to  work  with  time-independent 
methods;  however,  some  authors  [45-50]  discuss  the  approach  as  it  applies  to  time 
dependent  methods.  The  method  we  use  is  similar  to  the  “close-coupled  wave 
packet”  (CCWP)  method  [45-47],  because  the  wavefunction  is  represented  as  a  function 
of  R,  while  the  angular  variables  are  replaced  by  their  conjugate  momenta,  described  by 
coupled  equations. 

A  difficulty  stems  from  the  consideration  of  angular  momenta  that  must  be  added 
together  to  obtain  the  total  angular  momentum.  These  include  the  rotational  angular 
momentum  of  the  H2  diatom,  the  orbital  angular  momentum  of  the  B  -  •  -  H2  pair,  and  the 
open-shell  electronic  angular  momentum  of  the  boron.  Other  angular  momenta  that  are 
neglected  are  electronic  angular  momentum  of  H2  (it  is  assumed  to  be  in  a  Z  state)  and 
nuclear  spin  of  the  boron.  There  are  numerous  angular  momentum  coupling  schemes 
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(Hund’s  cases)  which,  while  equivalent  if  treated  fully,  are  intended  to  provide  alternative 
approximate  approaches  [33,51]. 


Space-Fixed  Coordinate  Systems 

We  are  primarily  concerned  with  identifying  the  location  of  the  nuclei.  Since 
there  are  three  atoms,  there  are  nine  degrees  of  freedom,  requiring  nine  coordinates.  In 
general  treatments,  the  three  atoms  are  typically  labeled  A,  B,  and  C,  where  for  our 
problem  A  labels  the  boron,  and  BC  labels  the  hydrogen  molecule.  The  first  coordinate 
system  that  may  be  considered  is  a  space  fixed  laboratory  coordinate  system,  shown  in 
Figure  1,  in  which  the  nine  coordinates  are  the  three  Cartesian  coordinates  of  the  three 
atoms.  The  Hamiltonian  can  be  transformed  to  a  space  fixed  center  of  mass  (SF) 
coordinate  system,  shown  in  Figure  2,  by  separating  out  the  coordinates  RCM  = 

( Xcm,  Ycm,  Zc, „ )  describing  the  center  of  mass  motion  of  the  system.  The  other  six 

coordinates  are  R  ,  the  vector  between  the  A  atom  and  the  center  of  mass  of  BC,  and  r  , 
the  vector  between  the  B  and  C  atoms. 

It  is  desirable  to  use  a  set  of  internal  coordinates  to  describe  the  system.  The 
potential  energy  surfaces  are  found  by  solving  the  Schrodinger  equation  in  a  system  in 
which  the  electronic  coordinates  are  measured  relative  to  the  nuclear  coordinates  -  a 
body-fixed  (BF)  coordinate  system.  The  internal  coordinates  that  are  used  are  often 
called  Jacobi  scattering  coordinates.  Let  A  label  the  boron,  and  BC  label  the  hydrogen 
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molecule,  as  is  typical.  Then  the  (A,  BC)  Jacobi  coordinates  are  R,  the  distance  from  the 
center  of  mass  of  the  BC  diatom  to  the  A  atom;  r,  the  distance  between  B  and  C,  and  9, 
the  angle  between  the  H2  bond  and  the  line  connecting  the  B  atom  to  the  H2  center  of 
mass.  One  can  also  use  the  (B,  AC)  or  (C,  AB)  Jacobi  coordinates;  the  three  possible 
sets  of  Jacobi  coordinates  label  three  possible  arrangement  channels.  For  our  problem, 
which  is  nonreactive,  the  system  stays  in  the  (A,  BC)  channel,  making  those  coordinates 
the  most  convenient.  If  we  use  the  (A,  BC)  internal  coordinates,  together  with  the  three 
coordinates  that  describe  the  center  of  mass  motion  of  the  system,  we  still  need  three 
more  coordinates  to  describe  the  orientation  of  the  ABC  system  in  space.  These  are 
found  by  transforming  to  a  body-fixed  (BF)  coordinate  system.  There  are  two  common 
ways  to  do  this: 


Figure  1 .  Space- fixed  coordinates 
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Figure  2.  Space-fixed  center  of  mass  coordinates 

Bodv-Fixed  Coordinate  Systems 

The  space  fixed  center  of  mass  coordinate  system  can  be  transformed  to  a  quasi¬ 
body  fixed  coordinate  system  defined  as  follows:  Locate  the  body-fixed  z-axis  to  point 
from  the  center  of  mass  of  the  H2  molecule  to  the  B  atom,  then  locate  the  body-fixed 
y-axis  in  the  space- fixed  xy-plane.  The  transformation  from  SF  to  BF  coordinates  can  be 
described  by  Euler  angles  (aft  0 ) ,  where  a  and  ft  are  the  azimuthal  and  polar  angles  of 
the  B  +  H2  rotor  with  respect  to  the  SF  frame.  The  internal  coordinate  6 is  the  polar  angle 
of  the  H2  diatom  with  respect  to  the  BF  frame,  and  (f  is  the  azimuthal  coordinate  of  the  H2 
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with  respect  to  the  BF  frame.  The  complete  set  of  body- fixed  coordinates  is  thus  Xcm, 
Ycm,  ZCm,  a,  R,  r,  6. '  and  (f> .  These  coordinates  are  illustrated  in  Figure  3,  and  will  be 
used  later  to  define  a  representation  for  the  Hamiltonian. 


x 

Figure  3.  Body- fixed  coordinate  system  I 


For  completeness,  another  body- fixed  coordinate  system  can  be  defined  in  a 
similar  fashion.  Locate  the  BF  z-axis  in  the  same  fashion  as  above.  Then  locate  the  BF  x 
axis  in  the  plane  of  the  ABC  system,  perpendicular  to  the  BF  z-axis.  The  transformation 
in  this  case  is  described  by  Euler  angles  (  a  f3  y) .  The  complete  set  of  body-fixed 
coordinates  is  thus  Xcm,  Ycm,  Zcm,  a,  ft,  y,  R,  r,  and  0,  and  is  illustrated  in  Figure  4. 
However,  these  coordinates  will  not  be  considered  further  in  this  work. 
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r 


In  Eq.  ( 10 )  pR  and  pr  are  the  momentum  operators  for  motion  of  A  relative  to  BC  and 
the  vibrational  motion  of  BC.  j  is  the  orbital  angular  momentum  operator  (BF  or  SF)  for 
the  BC  molecule,  and  L  is  the  orbital  angular  momentum  operator  for  A  and  the  BC  pair. 
Vel  and  Vso  are  the  electrostatic  and  spin-orbit  potential  operators. 


Restriction  to  Two  Internal  Degrees  of  Freedom 

To  restrict  the  problem  to  two  internal  degrees  of  freedom,  the  BC  bond  length  r 
will  be  held  fixed  at  req  .  Thus  the  p]  term  is  zero,  and  the  potential  will  be  evaluated  at 
the  fixed  length  req,  which,  for  H2,  is  taken  to  be  1 .402  atomic  units,  as  computed  by 
Alexander  [14].  This  approximation  amounts  to  neglecting  the  vibrational  degree  of 
freedom  of  the  hydrogen  molecule,  which  can  be  a  good  approximation  at  low  energies, 
since  the  threshold  for  the  v  =  0  -»  v-  1  vibrational  transition  for  H2  occurs  at  an  energy 
of  0.0199  atomic  units,  which  is  larger  than  the  highest  energies  for  which  we  obtain 
S-matrix  elements  (see  Section  V).  The  hydrogen  molecule  is  further  assumed  to  be  in  a 
'f  state,  so  that  it  contributes  zero  electronic  angular  momentum.  Table  1  lists  some  of 
the  constant  parameters  in  the  Hamiltonian. 
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Table  1 .  Parameters  in  the  Hamiltonian 


Name  Meaning  Value  (atomic  units') 


mH 

Mass  of  hydrogen 

1837.15 

mB 

Mass  of  boron  1 1 

20208.7 

Mbc 

Reduced  mass  of  H2 

918.58 

Ma.bc 

Reduced  mass  of  B  -  •  -  H2  pair 

3109.03 

€ 

Spin-orbit  constant  of  boron 

4.876  x  10 

V  eo 

Fixed  H2  bond  length 

1.402 

Choice  of  Basis 


Once  the  Hamiltonian  is  identified  in  operator  form,  it  must  be  represented  in  a 
basis  to  be  considered  numerically.  A  good  basis  is  identified  by  Dubemet  and  Hutson  in 
a  paper  that  outlines  several  possibilities  [33].  We  use  their  “Case  1A”  angular  basis  set, 
which  is  named  for  Hund’s  coupling  case  “A”.  The  Case  1 A  coupling  diagram  is  shown 
in  Figure  5  (from  Ref.  [33]).  This  coupling  scheme  in  which  ja  ,j  and  /project  separately 
onto  the  BF  z-axis  provides  a  simple  framework  to  neglect  (if  numerically  justified)  some 
of  the  coupling  terms  between  these  angular  momenta  which  arise  from  the  angular 
kinetic  energy  terms  in  the  Hamiltonian,  but  not  the  potential  energy.  This  is  done  later 
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in  Section  V,  through  the  centrifugal  sudden  approximation.  With  some  notational 
changes,  the  full  basis  is 


JJlsjaR)  (12) 

Mk  co  / 

In  this  basis  J is  the  total  angular  momentum  quantum  number  and  M its  projection  on  the 
space-fixed  z-axis.  j  and  k  are  the  diatom  orbital  angular  momentum  and  its  projection 
on  the  body-fixed  z-axis.  /  and  5  are  the  orbital  and  spin  angular  momenta  of  the  boron 
atom  and  take  the  values  1  and  1/2  for  the  2 p  orbital.  ja  and  co  are  the  boron  atom  total 
angular  momentum  and  projection  on  the  body-fixed  z-axis.  R  is  the  Jacobi  coordinate 
described  above,  and  r  is  omitted  due  to  the  constraint  to  two  degrees  of  freedom.  Since  j 
and  ja  are  present  in  the  basis,  terms 


Figure  5.  Coupling  case  “1  A” 
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Table  2.  Angular  Momentum  Quantum  Numbers 


Quantum  Number 


Meaning 


J 

M 

P 

j 

k 

L 

l 

s 


ja 

CO 


Total  angular  momentum  (conserved) 

Projection  of  Jon  SF  z-axis  (conserved) 

Projection  of  Jon  BF  z-axis;  equal  to  k  +  co 

Diatom  (H2)  angular  momentum 

Projection  of  j  on  BF  z-axis 

Orbital  angular  momentum  of  B  +  H2  pair 

Electronic  orbital  angular  momentum  of  B  (unpaired) 

Spin  angular  momentum  of  B  (unpaired) 

Total  B  atom  electronic  angular  momentum 
Projection  of  ja  on  BF  z-axis 


Using  the  coordinate  system  where  the  BF  x-axis  is  in  the  SF  xy-plane  (BF  coordinate 
system  I),  this  basis  represented  in  angular  coordinates  is  [33] 


ap  d<j)\lsJa  R' 
co 


r2J  +  \_ 
y  Aft 


v  1/2 


J  J  lsja  R) 

M  k  co  I 

(a,p,0)  Yjk(0,p)  <5{R-R') 


(13) 


<nJ * 

■uM(k+a>) 
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where  k  +  co is  the  projection  of  J onto  the  BF  z-axis  (Z  can  have  no  projection  on  the  BF 
z-axis),  and  is  called  P  in  Ref.  [33].  The  rotation  matrix  element  multiplied  by  the 
prefactor  is  the  rotor  wave  function  for  the  A-BC  system  that  describes  its  orientation 
with  respect  to  the  SF  axes  with  the  angles  a  and  (3.  The  spherical  harmonic  is  the 
wavefunction  for  the  BC  diatom,  which  describes  its  orientation  with  respect  to  the  BF 
axes,  with  angles  <9  and  cj). 


Representation  of  Each  Term  of  the  Hamiltonian 

Using  the  above  basis,  we  evaluate  matrix  elements  of  the  Hamiltonian.  For 
notational  convenience,  we  omit  J ,  M,  l,  s,  and  R  from  the  basis  vectors,  since  J  and  M 
are  conserved  by  the  Hamiltonian;  also  we  only  consider  2Pj  states  of  boron,  so  /  and  s  are 
conserved  by  our  Hamiltonian  (though  this  is  only  approximate).  R  is  implicitly  assumed 
to  be  present.  Matrix  elements  of  the  Hamiltonian  are 


fjjJa  <» 
n  jf  kk' 


k  co 


f  Ja\ 

k' '  CD') 


(14) 


Matrix  elements  for  each  term  in  the  Hamiltonian  can  be  considered  separately.  For  the 
radial  kinetic  energy  term,  the  matrix  elements  are 
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IJJo  Pi 

\k  ®  2 Ma,bc 


— — — - — ~~jR  8* 
2 HAfiC  R  dR2 


(15) 


8 xx'  -  8jj'MM'IVss'  8 jj'kk'  8]ja>m»'  8{R  R  ) 


where  the  Sxx,  shorthand  indicates  that  the  term  is  diagonal  in  all  the  angular  quantum 

numbers  (Though  a  delta  function  in  R  is  present,  the  kinetic  energy  is  of  course  not 
diagonal  in  R,  due  to  the  derivative).  The  angular  kinetic  energy  of  the  BC  diatom  is  also 
diagonal: 


2  Hi 


r2 

BC  req 


n  2  ^XX' 

2PBC  req 


bj(j  +  l)8xx. 


(16) 


where  the  above  equation  also  defines  the  rotational  constant  b  of  the  H2  rotor.  However, 
the  L2  term,  the  “tumbling”  energy  of  rotation  of  the  BC  diatom  and  the  A  atom 
together,  is  not  diagonal  and  poses  some  difficulty.  This  term  must  be  evaluated  by  using 
the  fact  that  L2  =  (J  -  j -  ja)2  : 

Aa  a  .  A*  A*  A  A  A  A  A  A 

L  =J  +  j  +ja  -  2J  ■  j  -  2J  •  ja  +2 j'ja 

-A  „  A«  /A/,  A  A  A  A  A  A 

=  J  +j  +Ja  -J  j-~J  j+~2Jzjz 

A  A  A  A  A  A  11/  ) 

-J+ja.-J-ja+-2JJaz 

A  A  A  A  A  A 

+  hi  a-  +  JJa+  +2  jJaz 
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The  diagonal  portion  is  thus 


I  JJa 

\k  CO 

-  - . -z  ( /(/  + 1) ■ +  j(j  + 1)  +  ja Ua  +  l)-2(k  +  o)f+2ko)) 5* 

2 Ma,bc  R 


^Ma.bc  R 


while  there  are  off-diagonal  terms  in  k,  co,  and  P  =  k+  co,  arising  from  the  raising  and 
lowering  operators.  The  actions  of  the  terms  in  Eq.  ( 17  )  that  involve  raising  and 
lowering  operators  on  the  basis  functions  are 


A  ,  A 

H 


aT 


7  ±J a 


77a' 
k  coj 

JJa ' 

k  0)j 

77  a' 
k  G)i 


=  n2(j(J  +  l)-P(P  +  l))m  (7(7  + 1) _ k(k  + 1))1/: 
=  n2  {J(  J  + 1)  -  P(P  + 1)),/2  0(7  + 1)  -  0)(0)  + 1))1 
=  (7(7  + 1)  -  i  I))*  2  (7(7  + 1)  _  + 1))’  2 


7  7a\ 

A:  +  l  0)  j 

j  Ja  \ 
k  co  +  H 


(19) 


7  7a  ' 
k±  1  <27  +  1/ 


The  raising  and  lowering  operators  j±  and  ja±  behave  in  the  usual  manner,  but  the 

operators  J±  behave  differently,  since  the  body-fixed  operators  must  be  used  fori , 
denoted  by  the  raised  ±  signs.  The  actions  of  the  body- fixed  raising  and  lowering 

operators  are  interchanged  [55]  so  thati+  is  a  lowering  operator  andi"  is  a  raising 
operator,  and  both  affect  the  quantum  number  P  of  the  state.  Hence  in  the  three  terms 
considered  above,  the  action  of  the  two  operators  is  first  to  apply  a  raising  (or  lowering) 
operator  to  raise  (lower)  k  or  co.  In  the  first  two,  the  body-fixed  raising  (lowering) 
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operators  then  raise  (lower)  P  to  match  the  fact  that  k  or  co  has  been  raised  (lowered)  -  for 
a  ket  to  be  valid,  P  must  equal  k  +  co .  In  the  last  case  in  Eq.  ( 19  ),  the  other  of  k  or  co  is 
lowered  (raised). 

Thus  the  kinetic  energy  term,  through  the  contribution  of  the  L2  term,  is  not 
diagonal  in  the  total  body-fixed  projection  P  =  k+  co.  However,  the  overall  magnitude  of 
the  L2  term  is  usually  small  for  our  problem,  for  low  values  of/,  due  to  the  small 
magnitude  of  the  rotational  “constant”  B  (R)  =  ft2  /(  2  juABC  R2 )  relative  to  the  potential 

energy  (see  section  V).  Later  the  terms  off-diagonal  in  P  will  be  neglected,  but  the  third 
portion  of  Eq.  ( 19 ),  which  is  diagonal  in  P  but  not  in  k  and  co,  will  be  kept. 

Next  we  require  the  matrix  elements  of  the  potential  in  this  basis.  First  consider 
the  electrostatic  potential  Vel .  Dubemet  and  Hutson  expand  the  potential  as  follows: 

vj&MiJ.)  =  Z  h  h  ( 2o ) 

where  Qm  is  a  renormalized  spherical  harmonic  [56],  6a  and  <j)a  are  electronic  coordinates 
of  the  boron  2p  orbital  with  respect  to  the  BF  z-axis,  and  the  sum  over  //ranges  from 
-  min  (  Ar ,  Aa  )  to  +  min  (Ar,Aa).  The  coefficients  (R)  are  obtained  from  the 

calculation  of  the  potential  energy  surfaces.  Since  the  potential  energy  for  the  electrons  is 
calculated  in  the  body  frame,  the  potential  is  actually  a  function  of  <f>-  <pa  rather  than  cj) 
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and  <pa  separately.  This  is  the  reason  that  the  same  //  appears  in  the  two  spherical 
harmonics  [33].  For  the  interaction  of  a  closed-shell  molecule  with  an  atom  in  a  P  state, 
only  terms  with  Aa  =  0  or  2  are  present  [33]. 

The  matrix  elements  of  the  electrostatic  potential,  in  the  above  expansion,  are 
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\k  CO 
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k'co'i 


(21) 


XrXapi 


c^(0,p) 


k'l 


Ja 

\CO 


fy 

Ja 
OJ  , 


JJ'MM 

ll'ss 


where 


\  =  (-1)*  [(2y  + 1)(2/+1)] 


1/2 


f  J  A  fsf 

-k  fj.  k' 


j  A  J 
0  0  0 


(22) 


and 


Ua 

\oo 


c^{0a,<Pa) 


Ja 


CO' 


,  =  (-irs(2/  +  l)[(2ya+l)(2y;+l)] 


1/2 


[Ja 

A  Ja  1 

'  j'a 

A 

J  a  " 

7 

K  0 

V 

5  /  J 

Kco' 

0  0, 

(23) 


29 


were  computed  using  the  Wigner-Eckhart  Theorem  and  Racah  algebra  formalism  by 
Dubemet  and  Hutson  [33].  The  symbols  in  large  parentheses  (:::)  and  braces  {:::}  are  3-j 
and  6-j  symbols  [56].  Since  a  3-j  symbol  is  zero  unless  the  lower  row  sums  to  zero,  the 
first  3-j  symbols  in  Eqs.-(  22 )  and  ( 23  )  require 


which  implies 


-k  +  jU  +  k'  =  0 
co'  -  ju-  co  =  0 


k  +  co  =  k'  +  a> 


(24) 


(25) 


This  reveals  that  the  electrostatic  potential  is  diagonal  in  the  total  body-frame  projection 
quantum  number  P  =  k+  co.  Additionally,  the  expansion  coefficients,  VXrXaM  (R)  , 

discussed  later,  will  be  nonzero  only  for  even  values  of  Ar .  This  is  a  consequence  of 
symmetry  arising  from  the  fact  that  H2  is  a  homonuclear  molecule  [14].  The  3-j  symbols 
have  the  property 


''a  b  c ' 

v°  0  °z 


(26) 


if  a  +  b  +  c  is  odd  [56].  Together  with  the  fact  only  even  values  of  Ar  are  present,  this 
means  that  the  second  3-j  symbol  in  Eq.  (  22  )  implies  that  the  transitions  will  take  place 
separately  within  the  sets  of  even  and  odd  diatom  rotational  quantum  numbers  j  . 
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The  final  term  in  the  Hamiltonian  is  the  spin-orbit  term  Vh .  The  usual  form  of  this  term, 

/( R )  l-s ,  requires  the  dependence  of  the  spin-orbit  coupling  on  the  distance  R.  It  is  a 
common  assumption  (the  “pure-precession”  approximation  [14]),  which  we  make,  to 
neglect  the  R  dependence  of  this  term,  assuming  that  the  interaction  does  not  perturb  the 
atomic  orbitals  significantly.  The  spin-orbit  term  is  diagonal  in  our  basis  with  matrix 
elements 


IJJa 
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(27) 


where  the  value  of  ^for  boron  is  10.7  cm'1,  or  4.876  x  10  atomic  units  [14,27]. 


Potential  Expansion  Coefficients 

One  still  needs  the  specific  form  of  the  potential  energy  expansion  coefficients 
u  C^)  to  compute  the  potential  matrix  elements.  These  coefficients  have  been 

determined  by  ab  initio  calculations.  Three  adiabatic  potential  surfaces  describing  the 
B  +  H2  interaction  have  been  computed  as  functions  of  R  and  6  by  M.  Alexander  [14] 
using  multireference,  configuration-interaction  calculations.  In  these  calculations,  the  H2 
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bond  length  was  held  fixed  at  1 .402  atomic  units  as  discussed  earlier.  The  surfaces 
correspond  to  the  three  degenerate  atomic  2p  orbitals,  where  spin  has  been  neglected. 
From  the  adiabatic  surfaces,  diabatic  potentials  Vxx  ,  Vyy,  Vzz  and  Vxz  have  also  been 
computed  by  Alexander  as  functions  of  R  and  fusing  an  approximate  transformation 
approach,  developed  by  Rebentrost  and  Lester  in  their  study  of  the  similar  F  +  FL 
interaction  [8,14],  The  Vxx ,  Vyy,  and  Vzz  diabatic  surfaces  nominally  correspond  to  the 
three  orientations  of  the  Cartesian  px  ,py ,  and  pz  orbitals  [14].  The  off-diagonal  coupling 
Vxz  corresponds  to  an  interaction  between  the  two  p  orbitals  in  the  triatomic  plane.  The 
2D  surfaces  may  be  obtained  from  HIBRIDON™  [57,58].  This  code  has  been  graciously 
provided  to  us  by  Alexander. 

In  Ref.  [14],  Alexander  expands  the  ^dependence  of  the  potential  energy  surfaces 
using  reduced  rotation  matrix  elements  d'm0  (ff) ,  where  /  and  m  are  arbitrary  indices. 

Since  the  third  index  is  zero,  these  are  equivalent  to  spherical  harmonics  [56]: 


4o(*)  = 


2k 
2/  +  1 


1/2 


Ylm(W)e 


(28) 


By  defining  Vs  and  Vd  to  be  one-half  the  sum  and  difference  of  the  Vxx  and  Vyy  potential 
energy  surfaces,  so  that 


vxx  =  (K-vd) 

Vyy  =  (Vs+Vd) 


(29) 
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the  angular  expansion  was  made  by  Alexander  as  follows: 


V22(R,&)  =  y£vr(R)d‘0(^ 
1=0 

VAR.#) 

/=0 

VAR.0) 

1=2 

v^R,e)  =  Y/rw4 oW 


The  independent  expansion  coefficients  VZZ{R),  Vs (R)  ,V,d(R),  andF/“(i?)  are 
described  in  Ref.  [14]  and  are  available  in  the  HIBRIDON™  code.  For  an  interaction 
with  a  homonuclear  molecule  such  as  H2,  the  odd-/  expansion  coefficients  vanish  due  to 
symmetry  [14].  The  above  expansion  coefficients  are  connected  to  the  BF  expansion 
coefficients  (R)  that  we  need  for  Eq.  ( 21  )  by  the  formulas  given  in  Ref.  [15] : 


VXr00(R)  =  [2Vl(R)  +  V^(R)]/3 
^20  (R)  =  5[-Vl(R)  +  vZ(R)]/3 

v^(R)  =  srz(Ry& 


Note  that  VKX  M  (R)  =  VA^ (R)  [33].  Since  the  odd-/  coefficients  in  Eq.  (  30  )  are  zero, 

the  BF  expansion  coefficients  in  Eq.  (  3 1 )  are  thus  zero  for  odd  values  of  Ar ;  this 
satisfies  the  requirement  for  the  conclusion  earlier  that  the  H2  rotational  transitions  will 
take  place  separately  within  the  sets  of  odd  and  even  j  quantum  numbers. 
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Since  J  and  M  are  conserved,  we  first  choose  a  block  of  the  Hamiltonian  matrix  labeled 
by  J  and  M .  An  approximation  must  be  made  later  to  truncate  the  (infinite)  matrix  by 
only  including  basis  vectors  up  to  a  maximum  value  ymax  of  the  diatom  rotational  quantum 
number  j  .  This  may  or  may  not  be  justifiable  by  the  energetic  accessibility  of  the  higher 
j  rotational  states  of  the  diatom,  the  energies  of  which  are  listed  in  Table  3.  If  the  wave 
packet  energies  are  lower  than  the  cutoff  needed  to  make  a  transition  to  a  higher  j,  energy 
conservation  forbids  that  state  from  being  accessible. 


Table  3.  Rotational  levels  of  H2 


Energy  (atomic  units) 


0 

1 

2 

3 

4 

5 

6 
j 


0 

0.00055 

0.00166 

0.00332 

0.00554 

0.00831 

0.01163 

2.769  x  10'4  x  j  0+1) 
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IV.  Implementation 


Numerical  Representation  of  Wave  Packets 

To  numerically  represent  a  wavefunction,  so  that  it  can  be  propagated  in  the 
application  of  the  channel  packet  method,  we  use  a  uniform  spatial  grid  of  discrete  points 
at  which  the  value  of  the  wavefunction  is  tabulated.  The  spatial  grid  must  be  large 
enough  in  extent  to  support  the  entire  wavefunction,  which  is  possible  because  the 
wavefunction  is  chosen  to  be  a  wave  packet  finite  in  extent.  The  discretization  of  a 
wavefunction  on  a  grid  of  Appoints  is  thus  given  by  N  numbers: 


=  *(*)  = 


'  *,(*,) ' 
^2(x2) 


(32) 


In  a  radial  or  Cartesian  direction,  the  Fourier  transform  grid  method  that  is  used  imposes 
requirements  on  the  grid.  The  periodic  boundary  conditions  assumed  by  the  Fourier 
transform  must  be  satisfied  by  the  wave  packet  approaching  zero  at  the  grid  boundaries, 
which  must  be  monitored  during  the  calculation.  The  same  must  be  true  of  the 
momentum  representation,  which  imposes  a  requirement  on  the  grid  spacing  in  the 
coordinate  representation.  The  maximum  momentum  representable  on  a  coordinate  grid 
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with  spacing  Ax  is  given  by  the  Nyquist  value  1/(2  Ax).  Additionally,  a  true  fast  Fourier 
transform,  scaling  as  O  (N  log  N)  instead  of  O  (A/2)  is  only  available  if  the  number  of  grid 
points  A  is  a  power  of  two;  if  more  coordinate  space  is  required  to  support  the 
wavefunction,  the  number  of  grid  points  must  be  doubled,  or  the  spacing  Ax  increased 
(which  reduces  the  maximum  representable  momentum). 


Gaussian  Wave  Packets 


A  one-dimensional  Gaussian  wave  packet,  at  t  =  0  with  initial  average  position  x0, 
initial  average  momentum  ko,  and  initial  position  uncertainty  Ax  is  defined  in 
Merzbacher  [59]: 


¥(*,  t  =  0)  = 
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(2kAx  ) 
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4  Ax2  0 
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To  analytically  evolve  this  wave  packet  to  time  t  =  r,  we  use  the  formula 


vp(x,  t  =  r)  = 
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The  initial  states  used  in  the  channel  packet  method  are  defined  to  be  Gaussian  wave 
packets.  This  allows  Eq.  (  34 )  to  be  used  to  analytically  propagate  the  wave  packet 
where  the  potential  is  zero  or  a  constant  value  (i.e.  during  the  creation  of  a  Moller  state). 


Propagation  of  the  Wave  Function 


The  key  element  in  a  time  dependent  technique  is  a  method  of  propagating  the 
wave  function.  Several  schemes  exist  and  have  been  reviewed  in  the  literature  [60];  we 
use  the  split  operator  method  [60].  This  technique  is  based  on  the  split  operator 

approximation  to  the  time  evolution  operator  U  (t)  =  exp  (-iHt/h) .  If  the  Hamiltonian 
for  the  nuclei  is  H  =  T+V ,  where  f  is  the  kinetic  energy  part,  and  V  is  the  potential 
energy  part,  one  can  write 


l^  +  AO)  = 


,  iV At.  iT  At  iV At  3 

exp( — — )  exp( - — )  exp( — — )  +  0(At  ) 


2  h 


h 


2  h 


l'P(O)  (35) 


The  term  0(Atf  arises  from  nonzero  commutators  involving  f  and  V ,  and  is  neglected 
by  using  a  sufficiently  small  time  step  At .  By  splitting  the  operator  into  separate 
exponentials  of  the  kinetic  and  potential  energy,  these  parts  can  be  evaluated  in 
representations  in  which  they  are  diagonal.  The  potential  energy  V ,  and  thus  its 
exponential,  is  diagonal  in  the  coordinate  representation  and  the  adiabatic  basis.  The 
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kinetic  energy  f ,  and  its  exponential,  are  diagonal  in  the  momentum  representation  and 
the  diabatic  basis. 

A  Fourier  transform  connects  the  coordinate  and  momentum  representations,  and 
a  unitary  matrix  transformation  connects  the  adiabatic  and  diabatic  bases.  An  application 
of  this  method  has  been  demonstrated  by  Alvarellos  and  Metiu  [61].  We  have  used  their 
example  wave  packet  propagation  to  test  the  code.  The  bulk  of  computational  time  is 
spent  performing  the  Fourier  and  matrix  transformations,  so  it  is  crucial  to  use  a  fast 
Fourier  transform,  which  requires  that  the  radial  grid  size  be  a  power  of  two.  The  unitary 
matrix  transformation  is  precomputed  by  one-time  diagonalization  of  the  potential  matrix 
at  all  values  of  R  .  Its  application  becomes  more  time-consuming,  and  it  requires 
increased  amounts  of  storage,  as  the  basis  size  increases.  The  application  of  the 
transform  requires  a  matrix-vector  multiplication,  but  not  a  diagonalization,  at  each  value 
of  R,  so  it  scales  both  in  time  and  storage  requirements  as  O  ( M2  x  N) ,  where  M  is  the 
basis  size  and  N  is  the  number  of  radial  grid  points. 


In  our  case,  the  radial  kinetic  energy  operator  becomes  -  — - in  the 

2^a,bc  R  dR 

coordinate  representation  and  it  is  most  convenient  to  store  the  “reduced”  wavefunction 

h2  d2 

R  lF(R) .  This  reduces  the  action  of  the  derivative  to - ,  which  is  evaluated 

2 jUAtBC  dR2 

by  transforming  to  the  momentum  representation  with  a  Fourier  transform.  When 
wavefunctions  are  integrated  to  obtain  the  correlation  function  (Eq.  (  6  )  in  the  channel 
packet  method),  we  must  remember  the  volume  element: 
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(36) 


C/r(t)  =  jyS'  (R)  exp (~iHt/h)^(R)  R2  dR 

0 

2 

Since  the  wavefunction  is  stored  as  R  5 F(R) ,  the  R  in  the  integration  is  already 
accounted  for.  Hence  the  same  computer  code  (representation  of  the  wave  packet  and 
numerical  integration)  can  be  used  for  spherical  coordinates  and  Cartesian  coordinates, 
but  the  interpretation  is  different. 


Absorbing  Boundary  Conditions 

Several  issues  can  arise  during  the  correlation  function  calculation.  The  grids 
used  to  calculate  the  Moller  states  may  be  different  from  the  one  used  to  calculate  the 
correlation  function.  An  example  of  this  stems  from  the  fact  that  absorbing  boundary 
conditions  (ABC)  can  be  used  to  significantly  speed  up  the  calculation  of  the  correlation 
function,  by  allowing  the  coordinate  grid  size  to  be  reduced.  If  the  evolving  Moller  state 
spreads  significantly  and  reaches  the  grid  boundary  before  it  has  ceased  to  overlap  the 
fixed  Moller  state,  it  will  reflect  off  the  edge  of  the  grid  due  to  the  periodic  boundary 
condition  nature  of  the  Fourier  transform  propagation  method,  unless  the  grid  is  made 
extremely  large  to  accomodate  the  full  wave  function.  To  avoid  this,  ABC  can  be  placed 
on  the  edge  of  the  grid.  The  ABC  take  the  form  of  an  imaginary  part  of  the  potential 
energy  added  to  each  adiabatic  (diagonal)  potential  surface: 
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(37) 


V  =  V,+iVABC 

where  V,  is  the  interaction  potential  and  VABC  is  the  absorbing  potential  [18].  When 
exponentiated  in  the  time  evolution  operator,  the  ABC  term  attenuates  the  wavefunction 
at  each  time  step  (if  VABC  has  the  correct  sign,  i.e.  negative  for  a  positive  time  step  At). 

Thus  if  ABC  that  do  not  overlap  the  fixed  Moller  state  are  placed  at  the  grid  boundaries, 
they  will  not  affect  the  value  of  the  correlation  function,  and  the  grid  size  can  be 
significantly  reduced.  For  a  good  example  of  the  application  of  ABC  to  a  two 
dimensional  problem,  see  Calfas  and  Weeks,  who  applied  it  to  the  test  case  of 
H  +  H2  [18,62].  The  absorbing  boundary  conditions  are  not  perfect;  they  always  cause 
some  reflection.  This  is  minimized  by  proper  choice  of  the  parameters  used  to  define 

VABC .  We  follow  [18]  and  use 


VABC  (R)  =  A  exp 


B 


(38) 


where  A  and  B  are  adjustable  parameters,  and  Rmax  is  the  value  of  R  at  the  grid  boundary. 
A  is  chosen  to  have  the  smallest  value  that  completely  attenuates  the  evolving  wave 
packet,  while  at  the  same  time  B  is  chosen  to  have  the  largest  value  that  keeps  the  ABC 
separate  from  the  initial  state  and  the  interaction  region.  These  choices  are  made  to  make 
the  ABC  as  smooth  and  gradual  as  possible,  to  minimize  reflection.  For  this  problem,  a 
suitable  choice  was  found  by  trial  and  error  to  be  A  =  0.0015  atomic  units,  and  B  =  30.0 
atomic  units. 
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Implementation  of  the  Channel  Packet  Method 


A  general  implementation  requires  three  separate  propagation  stages: 

I.  Preparation  of  the  reactant  Moller  state 

II.  Preparation  of  the  product  Moller  state 

III.  Evaluation  of  the  correlation  function 

In  stage  I  (stage  II),  the  initial  reactant  (product)  states  are  propagated  backward 
(forward)  in  time  under  the  asymptotic  Hamiltonian  to  a  time  t  sufficient  to  exit  the 
interaction  region.  They  are  then  numerically  propagated  forward  (backward)  in  time 
under  the  full  Hamiltonian,  back  to  t  =  0,  and  saved  to  disk  as  the  Moller  states. 

In  stage  III,  the  reactant  and  product  Moller  states  are  loaded  into  memory,  and  the 
reactant  state  is  propagated  backward  and  forward  in  time  until  the  correlation  function 
decays  to  zero. 

The  specific  problem  in  this  thesis  has  been  approached  using  explicit 
propagation  of  only  one  spatial  variable,  R  ,  since  the  problem  was  reduced  to  a  one¬ 
dimensional  propagation  on  many  surfaces  in  Section  III.  In  this  case,  we  can  justify 
(with  regard  to  computing  time)  a  simplification  of  the  above  procedure.  By  preparing 
the  Moller  states  in  the  asymptotic  limit,  the  Moller  states  are  equal  to  the  initial  states, 
and  stages  I  and  II  do  not  need  to  be  performed  (if  they  were,  the  resulting  Moller  states 
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would  be  identical  to  the  initial  states.)  Hence  only  stage  III  needs  to  be  performed,  and 
the  Moller  states  are  simply  defined  as  the  initial  reactant  and  product  wave  packets. 

Additionally,  in  the  computation  of  S-matrix  elements  for  transitions  from  the 
same  reactant  state  to  all  possible  product  states,  the  reactant  wave  packet  propagation  is 
identical.  Thus,  the  correlation  functions  between  the  evolving  reactant  Moller  state  and 
all  the  desired  product  Moller  states  are  computed  during  one  propagation  of  the  reactant 
Moller  state.  After  the  propagation,  each  correlation  function  is  used  in  Eq.  (  7  )  to 
compute  the  corresponding  S-matrix  elements. 


Accuracy  of  Wave  Packet  Propagation 

There  are  several  parameters  that  affect  the  numerical  accuracy  of  a  wave  packet 
propagation.  These  include  the  time  step  At,  the  coordinate  grid  spacing  Ax,  and  the 
coordinate  grid  size,  controlled  by  the  grid  spacing  and  the  number  of  grid  points  Nx.  The 
formal  error  in  Eq.  (  35  )  for  the  split-operator  method  is  O  (Hr).  The  numerical 
properties  of  the  split-operator  method  have  been  treated  [60].  A  key  positive  point  is 

2 

that  the  split  operator  method  is  unitary,  and  thus  norm-preserving  (the  integral  of  |  ¥  | 
stays  constant).  The  time  step  At  must  be  chosen  small  enough  so  that  the  product  of  the 
potential  or  kinetic  energy  and  At  is  also  small.  An  insufficiently  small  time  step  can 
result  in  dramatically  obvious  errors  in  the  wave  packet  propagation.  The  fundamental 
method,  however,  to  check  the  time  step,  is  to  compare  a  wave  packet  propagation  with 
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time  step  At  to  one  with  a  reduced  time  step,  say  At  12,  and  evaluate  the  convergence  of 
the  propagated  wavefunction  with  respect  to  the  reduced  time  step. 

The  nature  of  the  Fourier  transform  propagation  method  imposes  several 
conditions  on  the  spatial  grid  parameters  Nx  and  Ax .  First,  the  Fourier  transform  requires 
periodic  boundary  conditions.  If  a  wave  packet  approaches  one  end  of  the  grid,  it  will 
continue  by  wrapping  around  to  the  other  side,  if  there  is  no  potential  to  stop  it.  To 
satisfy  the  periodic  boundary  conditions  but  still  represent  the  true,  non-periodic,  physical 
system,  the  wavefunction  must  smoothly  approach  zero  at  the  grid  boundaries.  For  our 
problem  of  a  wave  packet  traveling  on  a  potential  surface  where  the  only  coordinate  is 
the  radial  coordinate  R,  the  potential  is  steeply  repulsive  at  the  small  R  edge  of  the  grid, 
and  the  wave  packet  is  defined  for  energies  that  cannot  (classically)  penetrate  the  barrier, 
so  it  effectively  decays  to  zero  on  that  edge.  However,  if  the  wave  packet  begins  to 
approach  the  large  R  grid  boundary,  it  will  incorrectly  reflect  from  the  boundary,  due  to 
the  steep  wall  at  the  small  R  edge  and  the  periodic  boundary  conditions.  Thus  during  the 
evaluation  of  the  correlation  function,  when  part  of  the  wave  packet  reaches  the  edge 
long  before  other  portions  have  exited  the  interaction  region,  absorbing  boundary 
conditions  are  used  to  prevent  the  wavefunction  from  reaching  the  grid  boundary. 
However,  the  absorbing  boundary  conditions  can  also  introduce  error  through  unintended 

reflections,  which  are  minimized  by  the  proper  choice  of  VABC  above. 

The  other  significant  issue  arising  from  the  nature  of  the  Fourier  transform  is  the 
analogous  requirement  that  the  momentum  grid  satisfy  similar  conditions  above.  The 
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momentum  grid  is  defined  by  the  choice  of  coordinate  grid.  It  has  a  maximum 
representable  momentum  kmax  =  n/  Ax,  (the  Nyquist  frequency)  which  is  entirely  defined 
by  the  coordinate  grid,  and  must  be  large  enough  in  extent  to  support  the  entire 
momentum  representation  of  the  wave  packet.  (In  atomic  units,  momentum  p  =  fik  =  k). 
If  the  wave  packet  has  momentum  components  larger  than  kmax,  those  components  will  be 
“aliased”  back  into  the  zone  with  |  k  \  <  kmax .  Thus,  during  the  propagation,  the 
momentum  representation  should  also  be  monitored  to  make  sure  that  it  does  not 
approach  the  momentum  grid  boundary. 

The  Fourier  transform  implementation  itself  is  unique  in  that  it  introduces 
extremely  small  amounts  of  numerical  error.  Due  to  the  fact  that  a  uniform  spatial  grid  is 
used,  the  sines  and  cosines  are  orthogonal  on  the  discrete  grid  [63:512],  and  thus  the 
discrete  Fourier  transform  (DFT)  performed  on  a  uniform  grid  is  an  exact  representation 
of  the  original  function  (on  the  uniform  grid).  The  operation  of  the  forward  and  inverse 
FFT  can  be  performed  millions  of  times  in  double  precision  and  only  change  the  very  last 
digits  stored  in  the  floating  point  representation,  due  to  computer  roundoff  error  only. 

To  compute  the  correlation  function,  the  overlap  integral  between  the  evolving 
reactant  Moller  state  and  the  product  Moller  state  is  computed,  using  the  trapezoid 
rule  [64],  which  has  an  accuracy  O  (Ac2) .  Since  the  integrand  approaches  zero  at  the 
endpoints,  the  trapezoid  rule  is  simply  evaluated  by  summing  the  values  of  the  integrand 
at  all  points  for  which  it  is  tabulated.  This  is  the  simplest  and  fastest  reasonable 
numerical  quadrature  technique  for  this  problem. 
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Accuracy  of  Computation  of  S-Matrix  Elements 


The  evolving  Moller  state  will  eventually  completely  exit  the  interaction  region, 
since  it  has  no  bound  state  components.  Once  it  leaves  the  interaction  region  and  also 
passes  the  region  of  the  fixed  Moller  state,  the  correlation  function  will  decay  to  zero.  At 
this  point,  the  evaluation  of  the  correlation  function  can  be  terminated.  Then,  the  Fourier 
transform  of  the  correlation  function  is  computed  and  used  in  Eq.  ( 7  )  to  obtain  the 
S-matrix  elements.  If  the  evaluation  of  the  correlation  function  is  terminated  before  it 
stabilizes  to  (numerical)  zero,  the  Fourier  transform  will  be  corrupted  with  oscillatory 
behavior  superimposed  on  the  true  answer.  This  can  be  seen  if  one  computes  a  “running” 
S-matrix  that  is  evaluated  over  time  as  the  correlation  function  is  computed. 

The  process  of  dividing  by  the  momentum  space  expansion  coefficients  7  (k)  is 
numerically  risky  for  values  of  k  where  7  is  small.  In  the  energy  ranges  where  the 
expansion  coefficients  are  small,  the  Fourier  transform  of  the  correlation  function  is 
usually  small,  but  contains  noise.  The  expansion  coefficients  decay  exponentially,  and 
thus  the  division,  and  the  S-matrix  elements,  blow  up  outside  the  range  over  which  the 
expansion  coefficients  are  significantly  nonzero.  Thus  the  answer  becomes  clearly  invalid 
outside  a  certain  energy  range,  since  the  absolute  value  squared  of  S-matrix  elements, 
which  represent  probabilities,  must  be  less  than  one.  This  limitation  in  the  valid  energy 
range  is  reasonable  -  we  should  only  get  useful  information  about  S-matrix  elements  for 
the  energy  range  that  we  include  in  our  wave  packets  that  explore  the  interaction 
potential.  This  problem  is  minimized  by  choosing  the  initial  wave  packets  to  have 
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expansion  coefficients  that  are  significant  in  the  entire  range  of  interest,  but  it  can  never 
completely  be  eliminated.  Another  related  problem  is  that  the  expansion  coefficients 
7  (k)  may  be  slightly  inaccurate  because  the  initial  wave  packets,  which  are  Gaussians  in 
coordinate  space,  are  never  truly  in  the  asymptotic  limit,  nor  are  they  entirely  separated 
from  the  absorbing  boundary  conditions.  This  can  cause  the  effective  wave  packet  that 
samples  the  potential  to  have  slightly  different  expansion  coefficients  than  those  used  in 
its  initial  definition.  This  problem  is  reduced  only  by  using  larger  and  larger  grids  to 
place  the  initial  state  closer  and  closer  to  the  “true”  asymptotic  limit. 


Validation  of  S-matrix  Elements 


S-matrix  elements  should  have  absolute  values  squared  that  are  less  than  one,  as 
discussed  earlier.  They  should  also  be  zero  at  energies  for  which  the  specific  transition  is 
energetically  impossible  (below  threshold).  One  way  to  test  the  numerical  parameters  in 
the  S-matrix  calculations  is  to  artificially  uncouple  the  potential  surfaces,  by  setting  the 
off-diagonal  terms  in  the  diabatic  potential  to  zero,  and  compute  S-matrix  elements  for 
reflection  where  the  reactants  and  products  are  on  the  same  surface.  The  resulting 
S-matrix  elements  should  be  exactly  one  for  all  energies.  The  energy  range  over  which 
this  test  numerical  answer  is  one  gives  an  indication  of  the  energy  range  over  which  the 
true  S-matrix  elements  should  be  valid,  which  should  correspond  to  the  numerically 
significant  range  of  the  momentum  expansion  coefficients.  In  Section  V  this  test  is  done 
for  the  parameters  used  in  our  calculation. 
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One  possible  qualitative  test  of  the  validity  of  the  S-matrix  elements  of  a  certain 
small  energy  range  is  to  propagate  a  wave  packet  defined  only  over  a  narrow  energy 
range  (so  that  it  is  very  broad  in  coordinate  space),  and  examine  the  “amount”  of  the 
wave  function  that  couples  to  the  product  state  in  question,  computed  by  integrating  the 
wave  packet.  This  value  is  then  compared  to  the  value  of  the  S-matrix  computed  by  the 
channel  packet  method  at  the  corresponding  energy.  This  qualitative  test  is  only  good  for 
S-matrix  elements  that  are  smooth,  broad  functions  of  energy,  when  one  wants  to  get  a 
general  idea  if  the  code  is  approximately  correct.  This  was  done  in  the  two  surface  case, 
as  an  initial  test  of  the  computer  code. 

A  strong  test  for  validity  is  the  fact  that  the  sum  over  all  possible  final  states  of 
the  probabilities  for  transition  from  a  given  initial  state  should  be  unity,  which  follows 
from  the  unitarity  of  the  S-matrix  [1 1].  Thus,  when  the  transition  probabilities  (absolute 
value  squared  of  S-matrix  elements)  are  computed  from  the  initial  state  to  all  the  other 
states  in  the  basis,  they  are  added  together,  and  the  range  over  which  this  sum  is  unity  is 
a  good  indicator  of  the  validity  of  the  S-matrices  over  that  range.  The  results  of  this  test 
will  be  presented  in  Section  V. 


Graphical  Visualization 

Code  to  implement  the  channel  packet  method  to  obtain  scattering  matrix 
elements  was  written  in  Fortran  90,  and  run  on  a  Silicon  Graphics  O2  computer.  An 
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important  part  of  the  implementation  is  the  ability  to  check  the  initial  conditions  and 
monitor  the  progress  of  the  calculation,  using  a  suitable  graphics  display.  Since  a  wave 
packet  propagation  may  take  hours  of  computer  time,  it  is  extremely  useful  to  be  able  to 
evaluate  the  initial  state  visually  to  ensure  it  is  correct  before  beginning  the  run.  Simple 
mistakes  such  as  preparing  the  wave  packets  in  the  wrong  location  (e.g.  off  the  grid)  or 
with  the  wrong  energy  components  become  immediately  obvious  upon  visual  inspection. 
As  the  wave  packets  evolve,  they  are  displayed  on  the  computer  screen,  together  with  the 
correlation  function  and  current  state  of  the  S-matrix  elements.  During  the  calculation, 
one  can  examine  these  to  ensure  that  the  calculation  is  still  running  as  planned,  or  to 
check  if  it  is  already  complete.  The  code  uses  OpenGL  for  its  graphics,  combined  with 
the  GL  Utilities  Toolkit  (GLUT).  OpenGL  is  a  platform  independent  3D  graphics 
application  programming  interface  based  on  the  Silicon  Graphics-specific  IRIS  GL. 

When  combined  with  GLUT,  which  provides  basic  windowing  ability,  it  allows  the  entire 
program  to  be  platform  independent,  and  not  tied  to  the  Silicon  Graphics  computer.  By 
use  of  modular  programming,  a  non-graphics  version  can  be  made  that  does  not  require 
the  OpenGL  library,  and  can  be  run  in  batch  mode. 

The  main  program  for  two  dimensional  calculations  is  controlled  by  a  data  file  of 
parameters  that  specify  exactly  how  the  initial  wave  packets  are  defined,  or  loaded  from 
files,  and  how  they  are  to  be  propagated  to  create  Moller  states  or  evaluate  the  correlation 
function.  The  potential  energy  surfaces  are  generated  by  a  Fortran  subroutine  provided  as 
part  of  the  HIBRIDON™  software  package  written  by  M.  Alexander,  who  graciously 
provided  us  with  the  subroutine.  HIBRIDON™  is  a  package  of  programs  for  the  time- 
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independent  quantum  treatment  of  inelastic  collisions  and  photodissociation  written  by 
M.H.  Alexander,  D.E.  Manolopoulos,  H.-J.  Werner,  and  B.  Follmeg,  with  contributions 
by  P.F.  Vohralik,  D.  Lemoine,  G.  Corey,  R.  Gordon,  B.  Johnson,  T.  Orlikowski,  A. 
Beming,  A.  Degli-Esposti,  C.  Rist,  P.  Dagdigian,  B.  Pouilly,  G.  van  der  Sanden,  M. 

Yang,  F.  de  Weerd,  and  S.  Gregurick. 

The  fundamental  graphics  subroutine  controls  the  plotting  of  a  function  of  two 
variables  as  a  grid  of  points,  lines,  or  a  smoothly  lit  surface,  which  can  be  scaled  and 
rotated  to  examine  it  from  all  angles.  This  plotting  subroutine  calls  commands  in  the 
OpenGL  API  (Application  Programming  Interface)  to  define  the  vertices  of  the  grid,  and 
the  lines  or  polygons  between  the  vertices  in  the  case  of  a  wire  mesh  or  smooth  surface. 
Separate  routines  control  color  and  lighting  parameters,  as  well  as  the  geometrical 
viewing  parameters  that  allow  the  function  plot  to  be  rotated  to  any  angle.  The  user  can 
select  any  function  (including  the  potential  surfaces,  wavefunctions,  correlation  function, 
and  S-matrix)  that  is  displayed,  and  through  a  keyboard  and  menu  interface,  the  user  can 
control  its  display  characteristics.  Each  plotted  variable  has  a  separate  set  of  display 
characteristics  stored  in  a  data  structure,  which  can  be  saved  to  disk  to  record  a  specific 
graphics  “state”  that  is  most  convenient  for  the  particular  application.  Through  the  use  of 
the  GLUT  API,  the  graphical  version  also  allows  interactive  control  over  the  process. 
GLUT  handles  the  window-system  dependent  tasks  such  as  opening  the  window  and 
receiving  input  via  mouse  clicks  and  pop-up  menus.  With  the  menu,  the  user  can  save 
and  load  wavefunctions,  correlation  functions,  and  S-matrices,  and  control  the  specifics 
of  the  evaluation  of  the  correlation  function.  After  the  calculation  is  complete  the 
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graphics  system  can  be  used  as  a  visualization  tool  to  examine  the  results.  For  intensive 
calculations,  the  graphics  can  be  used  first  to  check  the  initial  parameters,  and  then  turned 
off,  so  that  processing  time  is  not  allocated  to  unseen  graphics. 

In  the  calculations  where  the  basis  size  can  be  14  (or  larger),  the  radial 
dependence  of  the  wavefunction  is  a  14  element  column  vector  containing  14  “single¬ 
surface”  wavefunctions.  Rather  than  display  all  of  these  separately,  it  is  advantageous  to 
display  them  as  one  “two-dimensional”  wavefunction  where  one  coordinate  is  R  and  the 
other  coordinate  is  the  “basis  coordinate”.  With  this  approach,  the  entire  wavefunction 
can  be  examined  to  ensure  that  it  does  not  encounter  the  edge  of  the  radial  coordinate  or 
momentum  grids,  and  that  it  is  properly  absorbed  by  the  ABC.  One  can  also  examine 
how  much  coupling  has  taken  place  to  the  other  surfaces. 
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V.  Results  and  Discussion 


Truncation  of  Basis 


In  Section  III,  the  Hamiltonian  was  represented  in  matrix  form  using  a  basis  of  BF 
functions.  The  full  basis  is  infinite  in  size  and  must  be  truncated  to  treat  the  problem 
numerically.  Prior  to  truncating  the  basis,  one  can  take  advantage  of  the  fact  that  the 
Hamiltonian  is  block  diagonal.  Thus  one  can  select  specific  values  of  J  and  M  for  the 
problem,  and  treat  only  that  specific  block  (which  is  still  infinite  in  size). 

The  Hamiltonian  is  also  approximately  block  diagonal  in  the  total  body  frame 
projection  quantum  number  P  =  k+  co.  Off  diagonal  terms  in  P  arise  from  matrix 

elements  of  the  Z2  term.  An  approximation  that  we  will  make  is  to  neglect  these  terms. 
This  approximation  has  been  discussed  often  and  takes  many  names,  including  the 
centrifugal  sudden  (CS)  approximation,  helicity  decoupling  approximation,  or 
Jz-conserving  coupled  states  (also  CS)  approximation  [1 1,45-47].  If  the  magnitude  of  the 
rotational  constant  B(R )  of  the  A-BC  system  is  small,  due  to  a  large  reduced  mass 
combined  with  a  sufficiently  large  value  of  R  at  closest  approach,  and  the  total  J  is  small, 
this  may  be  a  good  approximation.  For  the  B  +  H2  problem,  the  distance  R  stays  at 
values  £  5  atomic  units,  so  the  magnitude  of  B  stays  small  relative  to  the  other  terms  in 
the  Hamiltonian. 
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Additionally,  as  stated  in  Section  III,  rotational  transitions  take  place  separately 
within  the  sets  of  odd  and  even  values  of  j,  the  rotational  quantum  number  for  the  H2 
molecule.  Hence  the  Hamiltonian  separates  into  two  blocks  labeled  by  the  sign  of  (-iy  , 
and  we  can  choose  one  to  work  in.  One  must  finally  truncate  the  basis  by  choosing  a 
value  jmax  that  is  the  highest  value  of  j  that  will  be  included  in  the  basis.  The  choice  of 
jmax  may  be  aided  by  energy  arguments;  however,  it  is  sometimes  necessary  to  include 
channels  of  higher  j  that  are  not  energetically  accessible  as  product  states  but  that  may  be 
temporarily  accessible  during  the  interaction.  For  completeness,  we  include  the  fact  that 
the  continuous  infinite  variable  R  is  replaced  by  its  finite  discretization,  discussed  in 
Section  IV. 


Pam  Basis  Sets 


For  this  problem,  then,  the  choice  must  be  made  for  the  values  of  J,  P,jmax ,  and 
whether  j  will  stay  even  or  odd.  A  straightforward  choice  for  an  initial  treatment  is  to 
choose  the  lowest  possible  value  of  J,  which  is  1/2,  and  choose  P  arbitrarily  to  be  +1/2  . 
The  Hamiltonian  matrix  is  independent  of  the  total  space-fixed  projection  M,  which  we 
can  leave  unspecified.  We  will  choose  to  study  transitions  starting  from  the  ground  state, 
SO  j  initial  will  be  zero,  and  hence  j  will  stay  even.  (The  H2  will  stay  para,  labeled  p-H2)  We 
study  several  choices  of  jmax,  and  hope  that  as  jmax  increases  beyond  the  energetically 
accessible  value,  the  scattering  results  converge  with  the  corresponding  increased  basis 
size. 
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Thus  our  basis  sets  will  be  labeled  by  the  available  values  of  j,  assuming  the  other 
choices  above.  The  simplest  basis  is  thus  denoted  by  (  0 ),  while  larger  bases  will  be 
labeled  ( 0,  2 ),  ( 0, 2,  4 ),  or  ( 0, 2, 4,  ...,jmax  ).  The  number  of  basis  vectors  in  this 
specific  basis  set  will  be  3  x  jmax  +  2.  The  reason  for  this  is  that  given  the  full  set  of  six 
spin-orbit  states  |  jaco)  of  boron,  each  will  have  a  specific  value  of  a>,  which  must  be 

paired  with  exactly  one  available  value  of  k to  give  the  value  P  =  k+  co=  1/2  .  This  is 
done  for  each  even  value  of y  except  for  j  =  0  ,  in  which  case  there  are  only  two  legal 
basis  vectors: 


j  J  a  \ 
k  coj 


(39) 


For  j  =  2,  there  are  six  basis  vectors: 


j  J'a\ 

k  coj 


€ 


2  K\  2K\ 

2  -%)’  l  -y2j 


2  A\  2  A\ 

1  -%/’  0  A/ 


(40) 


So  the  ( 0, 2 )  basis  contains  eight  basis  vectors,  combining  they  =  0  andy  =  2  basis 
vectors  listed  above.  There  are  again  six  basis  vectors  for  higher  values  of  eveny,  which 
can  be  seen  to  be  exactly  the  same  as  those  in  Eq.  ( 40 )  but  with  the  upper  left  “2” 
replaced  by  the  specific  value  of  j. 
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If  we  consider  matrix  elements  off-diagonal  in  P,  and  thus  include  all  possible 
values  of  P  in  the  basis  set,  it  will  be  significantly  larger.  Table  4  summarizes  the  size  of 
various  bases  containing  even  j,  which  we  denote  “ para ”,  with  and  without  making  the 
CS  approximation  (P  conservation),  with  J=  1/2. 


Table  4.  Size  of  para  bases  with  J=  1/2 


Basis  set 

CS  Basis  Size 

Non-CS  basis  size 

(0) 

2 

6 

(0,2) 

8 

36 

(0,  2,4) 

14 

90 

(0,2,4, 6) 

20 

168 

(  0,  2,  . . . ,  jmax  ) 

2  +  3  X  jmax 

6  x  X(2y  +  1) 

je  basis 


In  this  thesis  we  will  present  results  for  the  (  0 ) ,  ( 0,  2 ),  and  ( 0,  2,  4 )  bases, 
using  the  CS  approximation.  A  block  schematic  of  the  Hamiltonian  in  matrix  form  is 
shown  in  Figure  6.  The  darkest  shaded  area  represents  the  ( 0 )  basis;  this  together  with 
the  medium  shaded  area  represents  the  (  0,  2  )  basis  (8x8).  The  entire  shaded  grid 
represents  the  (  0,  2,  4  )  basis  (14x14 ).  Off-diagonal  terms  exist  throughout  the  entire 
matrix. 
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0  ja 
0  1/2 


2  ia 
k  co 


4  ja 

k  co 


Figure  6.  Block  schematic  of  the  Hamiltonian  matrix 


The  effective  potential  (potential  energy  +  angular  kinetic  energy)  matrix  in  the  ( 0, 2, 4 ) 
basis  has  142  =  196  elements  (though  it  is  symmetric).  To  present  the  potential  surfaces 
in  the  clearest  fashion,  Figure  7  shows  the  diagonal  elements  of  the  matrix  in  the  diabatic 
representation.  Figure  8  shows  the  first  row  of  the  potential  matrix  in  the  diabatic 
representation.  Figure  9  shows  the  adiabatic  potential  surfaces  ( the  diagonal  and  only 
non-zero  elements  in  the  adiabatic  representation). 
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R  ( atomic  units ) 


Figure  7.  Diagonal  elements  of  potential  matrix  in  diabatic  representation 


Figure  8.  First  row  of  potential  matrix  in  diabatic  representation 


Figure  9.  Adiabatic  potential  surfaces 


The  “One-Dimensional”  Problem 

The  even  j  rotational  levels  of  H2  label  />ara-hydrogen  (p-Eh),  with  nuclear  spin 
wavefunction  anti-symmetrized  and  spatial  wavefunction  symmetrized,  forming  a  singlet 
state.  The  odd  rotational  levels  are  ortho- hydrogen,  the  triplet  state.  At  room 
temperature,  hydrogen  is  found  in  a  ratio  of  approximately  3:1  ortho  to  para.  However, 
as  the  temperature  is  lowered  below  20  K  or  so,  the  hydrogen  molecules  freeze  into  their 
lowest  state  and  tend  to  stay  para,  with  frozen  hydrogen  being  mostly  para. 
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If  one  assumes  that  a  p-H^  molecule  in  its  ground  (j  =  0  )  rotational  state  stays  in 
the  ground  state  during  the  collision  with  a  boron,  the  ( 0 )  basis  set  in  Eq.  (  39 )  may  be 
used.  The  2x2  matrix  of  the  potential  Vel  +  Vso  in  this  basis  is 

h  =1/2  ja =3/2 

ja  =  1/2  (V0(R)+B(R)-f  -2 U2V2(R)  (41) 

ja=3/2{  -2  U2V2(R)  V0(R)+V2(R)+7B(R)+f/2; 

where  the  values  J=  1/2,  P  =  1/2  are  used  together  with  the  value  of  j„  in  Eq.  ( 18  )  to 
compute  the  coefficient  of  B(R) .  The  functions  Vo(R)  and  V2 (R)  are  given  in  Ref.  [14] 
from  combinations  in  the  expansions  in  Eq.  (  30  )  : 

V0{R)  =  [2VU{R)  +*£(*)]/ 3 

F2W  =  [-f,loW  +  V,T0(R)  ]  I  3 

From  Eq.  (  31 )  one  can  note  that  in  terms  of  the  body-fixed  potential  expansion 
coefficients,  V0(R)  =  Vooo(R)  and  5  F2(R)  =  V02o (R)  ■  The  radial  kinetic  energy  operator  is 
h2  1  d 2 

- -R  x  1  ,  where  1  is  a  2  x  2  unit  matrix.  Elements  of  the  2  x  2  matrix  in 

2Ma,bc  RdR 

Eq.  ( 41 )  are  shown  in  Figure  10,  together  with  the  adiabatic  surfaces  obtained  from  the 
matrix  diagonalization.  The  adiabatic  surfaces  are  the  largest  and  smallest  in  the 
interaction  region.  The  off-diagonal  coupling  term  approaches  negative  infinity  as  R 
decreases. 
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Figure  10.  Diabatic  (solid)  and  adiabatic  (dashed)  potential  surfaces  for  (  0 )  basis 

In  the  2  x  2  (j  =j  ’  =  0  )  case,  the  p-E^  can  be  considered  equivalent  to  a 
spherical  atom  colliding  with  the  boron,  and  the  mathematics  reduces  to  that  of  atom- 
atom  interactions,  which  are  somewhat  simpler.  The  three  basic  adiabatic  potential 
energy  surfaces  reduce  to  two,  due  to  the  extra  symmetry.  These  nominally  correspond 
to  the  p  orbital  of  the  B  atom  pointing  along,  or  perpendicular  to,  the  direction  of 
approach.  These  surfaces  are  one  dimensional  due  to  the  spherical  symmetry  and 
correspond  to  averages  of  the  isotropic  terms  in  the  full  potentials.  The  SF  basis  found  in 
Ref.  [33]  for  atom-diatom  interactions  becomes  simpler  in  the  atom-atom  case.  In 
Ref.  [14],  Alexander  uses  this  idea  to  evaluate  the  2x2  matrix  in  the  SF  frame  that 
results  when  j  is  constrained  to  stay  zero.  This  is  done  using  simpler  formulas  [65-66], 
The  resulting  matrix  is  the  same  as  we  find  in  the  body- fixed  case  except  for  a  sign  on  the 
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off  diagonal  terms,  and  slightly  different  forms  of  the  multipliers  in  front  of  the  B(R) 
terms.  This  difference  may  be  explained  by  our  use  of  the  CS  approximation. 
Numerically,  we  will  see  that  the  difference  in  the  form  of  the  angular  effective  potential 
has  little  effect  on  the  S-matrix  elements;  we  present  the  results  for  all  the  basis  sets 
using  the  BF  frame  for  consistency. 


Application  of  Channel  Packet  Method  to  “ID”  Problem 

To  compute  S-matrix  elements  for  the  B  (  2Pj /2 )  -»  B  ( 2P3/2 )  transition  in  the 
(  0  )  basis,  we  prepare  the  reactant  Moller  state,  using  Eq.  (  5  ),  from  an  initial  state 
( t  =  0 )  on  the  lower  (ja  =  1/2,  j  =  0,  labeled  2Pi/2  )  surface  and  the  product  Moller  state 
from  an  initial  state  on  the  upper  ( 2P3/2 )  surface.  The  initial  states  are  Gaussian  wave 
packets  entering  (reactants)  or  leaving  (products)  the  interaction  region.  Table  5  lists  the 
parameters  used  to  define  the  initial  states.  The  initial  momentum  and  momentum 
uncertainty  of  the  wave  packet  are  chosen  so  that  the  wave  packet  possesses  a  broad 
negative  (positive)  range  of  plane  wave  components  for  the  reactants  (products),  but  with 
no  appreciable  positive  momentum  (for  the  reactants)  or  negative  momentum  (for  the 
products).  As  stated  in  Section  II,  this  is  required  for  Eq.  ( 7  )  to  be  applicable.  The 
energy  range  is  also  chosen  so  that  no  plane  wave  components  penetrate  the  barrier  for 
small  R  in  the  potential. 
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Table  5.  Channel  Packet  Parameters 


Parameter 

Value  ( atomic  units) 

Reactant  initial  average  momentum 

-4.0 

Product  initial  average  momentum 

+4.0 

Reactant,  product  initial  position 

75.0 

Reactant,  product  initial  position  uncertainty 

0.65 

Correlation  function  time  range 

[  0 , 6.25  x  106  ] 

Coordinate  grid  spacing  Ax 

0.15 

Time  step  At 

25.0 

ABC  parameter  A 

0.0015 

ABC  parameter  B 

30.0 

The  first  frame  of  Figure  1 1  shows  the  reactant  Moller  state  at  t  =  0,  on  the  lower 
surface,  together  with  the  adiabatic  surfaces  and  absorbing  boundary  conditions  (dotted). 
The  remaining  frames  of  Figure  1 1  show  the  time  sequence  of  the  propagation  of  the 
reactant  Moller  state,  at  intervals  of  25,000  atomic  units.  The  wavefunction  on  the  lower 
surface  is  shown  with  a  solid  line,  while  the  wavefunction  on  the  upper  surface  is  shown 
with  a  dashed  line. 
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Figure  1 1 .  Propagation  of  reactant  Moller  state  in  (  0  )  basis 
from  t  =  0  to  t  =  300,000  atomic  units 


The  correlation  function  C3/2, 1/2  (0  in  Eq.  (  6  ) ,  shown  in  Figure  12,  is  computed 
by  holding  the  product  Moller  state  fixed  and  propagating  the  reactant  Moller  state  from 
t  =  0  forward  to  t  =  375,000  atomic  units.  The  numerical  limits  approximate  the  true  time 
interval  from  -00  to  +00,  and  are  chosen  to  include  the  range  over  which  the  correlation 
function  is  numerically  nonzero.  (In  the  ( 0 )  basis  the  correlation  function  has  already 
decayed  to  zero  by  this  time;  propagation  will  only  need  to  proceed  to  the  full  6.25  x  106 
atomic  units  when  the  larger  bases  are  used.)  When  the  reactant  Moller  state  reaches  the 
interaction  region,  some  of  it  begins  to  couple  onto  the  other  surface,  reflect  off  the 
barrier,  and  overlap  the  product  Moller  state,  resulting  in  a  nonzero  correlation.  The 
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absorbing  boundary  conditions  prevent  the  evolving  Moller  state  from  reflecting  from  the 
grid  boundary  after  it  has  passed  through  the  fixed  product  Moller  state.  The  value  of  the 
correlation  function  is  not  affected  by  the  attenuation  of  the  evolving  reactant  Moller 
state  in  the  area  of  the  ABC,  where  it  does  not  overlap  the  fixed  product  Moller  state. 


Figure  12.  Correlation  function  Czn,\n  it) 


The  Fourier  transform  of  the  correlation  function,  shown  in  Figure  13,  is  used  in  Eq.  (  7  ) 
to  compute  S-matrix  elements.  Also  required  are  the  momentum  space  expansion 
coefficients  T] ,  which  are  evaluated  implicitly  as  a  function  of  energy  using  Eq.  (  8  ),  and 
shown  in  Figure  14. 
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The  absolute  value  squared  |  S  f  gives  the  probability  PV2, 1/2  for  the  2Pi/2  — >  2P3/2 
transition,  and  is  shown  in  Figure  15.  The  transition  probability  exhibits  an  oscillatory 
dependence  on  energy.  This  is  similar  to  the  results  of  a  previous  study  of  transition 
probabilities  between  two  surfaces  using  a  close  coupled  method  [13],  in  which  there  is 
evidence  of  an  oscillatory  dependence  in  the  result  as  a  function  of  the  (discrete)  energies 
presented.  The  low  energy  portion  roughly  matches  the  shape  of  the  results  of  Flower 
and  Launay  for  fine  structure  transitions  in  collisions  of  C+  +  H2  [4],  who  constrained 
p-H2  to  the  j  =  0,  2  rotational  levels  in  the  calculation.  In  the  higher  energy  range,  higher 
rotational  levels  will  be  open  to  excitation,  which  requires  the  use  of  the  larger  (  0, 2 ) 
and  (  0,  2,  4 )  basis  sets. 


Figure  15.  Probability  P372’ 1/2  of  transition 
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To  compute  the  probability  for  no  transition  from  the  lower  state  ( 2P 1/2  — >  2Pi/2  ) 
the  product  Moller  state  is  prepared  on  the  lower  surface  rather  than  the  upper  surface. 
Over  the  numerically  valid  energy  range  where  the  initial  reactant  and  product  wave 
packets  have  nonzero  expansion  coefficients  Tj(  k  \E\ ),  the  sum  of  the  transition 
probabilities  for  the  reactions  starting  on  the  upper  surface  is  close  to  unity.  The 
significant  spikes  at  the  ends  mark  the  end  of  the  energy  range  for  which  the  expansion 
coefficients  are  significantly  nonzero.  There  is  a  small  overall  error  due  to  the  fact  that 
the  initial  states  can  never  be  truly  numerically  placed  in  the  asymptotic  limit,  which 
slightly  modifies  the  true  values  of  the  expansion  coefficients. 

The  threshold  energy  for  a  transition  from  the  lower  state  to  the  upper  state 
depends  on  the  energy  gap  between  the  states  and  the  behavior  of  the  coupling.  The  B 
spin-orbit  splitting  is  small  (7.31  x  10'5  atomic  units)  on  the  energy  scale  shown,  and  the 
numerical  results  are  in  question  in  this  region.  For  the  reverse  transition,  the  reactants 
are  prepared  on  the  upper  ( 2P3/2 )  surface  and  the  products  on  the  lower  ( 2P3/2 )  surface. 
The  probability  Pm’ 3/2  obtained  in  this  calculation  is  not  identical  to  Pil2, 1/2  ,  that  of  the 
reverse  reaction,  due  to  the  small  excitation  threshold. 

As  mentioned  in  Section  IV,  a  test  for  the  validity  of  the  S-matrix  element 
calculation  is  to  examine  the  case  of  “uncoupled”  surfaces  where  the  off-diagonal  terms 
in  the  diabatic  potential  are  set  to  zero.  In  this  case,  the  probability  obtained  from  the 
S-matrix  elements  for  reflection  from  a  given  surface  to  itself  should  be  unity,  over  the 
valid  energy  range,  (the  range  over  which  the  7  coefficients  are  numerically  significant). 
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Figure  16  shows  this  for  initial  states  on  the  lower  surface.  Here  the  valid  energy  range 
extends  up  to  energies  of  about  0.01  atomic  units. 


Figure  16.  Probability  of  reflection  pll2,m  in  the  uncoupled  case 


Application  to  Larger  Basis  Sets 

The  ( 0, 2 )  and  ( 0, 2, 4 )  basis  sets  in  the  CS  approximation  have  sizes  of  8  and 
14  basis  vectors  respectively  ( Table  4 ).  The  reactant  wave  packet  is  propagated  in  the 
same  manner  as  in  the  previous  discussion,  with  the  same  initial  parameters  (  Table  5  ), 
but  the  wavefunction  is  now  an  8  or  14  element  column  vector,  and  the  possibility  of 
coupling  to  the  higher  rotational  states  is  included.  The  correlation  functions  are 
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evaluated  between  the  evolving  reactant  Moller  state  on  the  lowest  surface  and  all  8  or  14 
possible  fixed  product  Moller  states. 

A  significant  difference  that  is  initially  observed  between  these  cases  and  the  ( 0 ) 
basis  is  that  the  correlation  functions  C3/2, 1/2  and  Cm,  1/2  computed  in  these  bases  are 
significantly  longer  lived  than  those  computed  in  the  (  0 )  basis,  and  require  evaluation  to 
6.25  x  106  atomic  units.  Shown  in  Figure  17  are  the  correlation  functions  C3/2, 1/2  where 
the  basis  is  ( 0,  2  )  and  (  0,  2, 4 )  (compare  to  Figure  12). 


74 


o 


0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 
0 

-0.01 

0  5  10s  1  106  1.5  106  2  10s  2.5  106  3  106  3.5  106  4  106 

Time  ( atomic  units ) 


(  0,  4 )  basis 


M/Wvww^^w- _ - 


J _ I _ I _ I _ I _ I - L 


Figure  17.  Correlation  function  Cm,  1/2  (t)  in  the  ( 0, 2 )  and  ( 0,  2, 4 )  bases 

The  long-lived  quality  of  these  correlation  functions  contributes  to  very  sharp  features  in 
their  Fourier  transforms,  and  thus  in  the  corresponding  S -matrix  elements,  shown  in 
Figure  18.  These  features  are  noticeable  primarily  in  the  energy  range  just  below  the 
threshold  for  they  =  0  -*j  =  2  rotational  transition.  The  reason  the  correlation  functions 
last  for  much  longer  is  that  wave  packet  components  that  have  almost  enough  energy  to 
make  the  0  -»  2  or  0  -» 4  transition  in  the  asymptotic  limit  will  actually  have  sufficient 
energy  in  the  interaction  region,  where  the  gap  is  smaller  due  to  the  potential  anisotropy. 
Hence  quasi-bound  states  can  be  formed  in  the  wells  of  the  rotationally  excited  potential 
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energy  surfaces.  These  states  are  never  fully  bound;  there  is  always  enough  energy  to 
escape  the  well,  but  it  may  take  a  very  long  time  to  couple  back  to  the  lower  surface. 
During  this  long  calculation,  it  is  essential  that  the  absorbing  boundary  conditions  cause 
absolutely  negligible  reflection  or  the  long  tail  of  the  correlation  function  will  be 
corrupted  by  contributions  from  reflections  by  ABC.  To  test  for  this,  the  calculation  was 
done  several  times,  with  different  values  for  the  ABC  parameters.  The  sharp  features 
remained  the  same;  additionally,  the  fact  that  the  sharp  features  add  to  unity  over  most  of 
their  range  provides  confirmation  that  ABC  reflections  are  negligible. 


(a)  (  0,  2  )  basis 
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(d)  ( 0,  2,  4  )  basis  magnified  to  show  sharp  features 
Figure  18.  Probabilities  of  transition  PV2’ 1/2  for  the  ( 0,  2  )  and  (  0,  2,  4  )  bases 

We  tabulate  the  transition  probabilities  for  the  reactants  starting  in  the  ground 
state  ( B  2Pi/2  and  j  -  0 )  to  the  available  product  states,  for  each  of  the  ( 0 ),  ( 0,  2 ),  and 
(  0,  2,  4 )  bases.  Thus  there  are  2,  8,  and  14  transition  probabilities  respectively.  In 
Figure  19,  we  show  the  sum  over  all  possible  final  state  transition  probabilities  starting 
from  the  given  initial  state  (the  ground  asymptotic  state).  Shown  in  Figure  20  are  the 
expansion  coefficients  //  of  the  initial  reactant  and  product  wave  packets  where  the 
reactants  are  prepared  in  the  ground  state  and  the  products  are  prepared  in  the  third  basis 
state  in  Eq.  ( 40  ):  (j  =  2,j„  =  3/2,  k  =  2,  6)=  -3/2 ).  The  deviations  from  unity  arise  in 
the  regions  in  which  at  least  one  of  the  rj  expansion  coefficients  that  contributes  to  an 


78 


S-matrix  element  in  the  sum  is  small,  and  from  the  fact  that  the  expansion  coefficients 
may  not  truly  represent  the  wave  packet  if  it  is  defined  too  close  to  the  interaction  region. 
This  latter  problem  can  only  be  reduced  by  using  larger  and  larger  grids.  The  transition 
probabilities  for  final  states  with  j  =  2  and  j  =  4  are  small  and  smaller,  respectively, 
indicating  that  the  results  are  converging  with  respect  to  increasing  basis  size. 

The  complete  set  of  transition  probabilities  from  the  ground  state  to  the  available 
product  states  for  the  three  bases  are  shown  in  Figure  21,  Figure  22,  and  Figure  23  in  the 
Appendix.  The  low  energy  portions  for  the  results  in  the  ( 0,  2  )  and  (  0,  2,  4 )  basis, 
below  the  j  =  0  -*j  =  2  transition  threshold,  are  similar  to  the  results  calculated  using  the 
( 0 )  basis  set,  but  the  sharp  features  are  different. 


(a)  ( 0 )  basis 
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Transition  probability  Transition  Probability 


(b)  (0,2)  basis 


(c)  (  0,  2,  4  )  basis 

Figure  19.  Sum  over  final  transition  probabilities  for  all  bases 
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Figure  20.  Expansion  coefficients  rj  of  initial  wave  packets  (product  has  j  =  2) 


Conclusion 


A  new  application  of  the  time  dependent  channel  packet  method  to  problems 
involving  coupled  electronic  surfaces  has  been  used  to  obtain  S-matrix  elements  for  the 
fine  structure  transition  that  occurs  in  collisions  of  B  and  H2 ,  including  the  possibility  of 
rotational  excitation  within  the  even  (para)  rotational  states.  The  CPM  gives  us  S-matrix 
elements  as  a  function  of  the  entire  energy  range  of  the  initial  channel  packets.  The 
probability  of  reaction  generally  has  an  oscillatory  dependence  on  energy,  and  for  low 
energy,  is  similar  in  form  to  that  of  a  previous  study  of  C+  +  H2  [4],  Sharp  features  show 


81 


up  in  the  S-matrix  elements  when  higher  rotational  states  are  considered,  due  to  the 
availability  of  temporary  potential  energy  for  rotational  excitation.  S-matrix  elements  for 
transitions  from  the  ground  state  to  all  final  states  in  bases  with  jmax  =  0,  2,  and  4  were 
tabulated. 

Recommendations  for  future  work  include  a  more  complete  study  of  S-matrix 
elements  given  the  framework  above,  such  as  considering  transitions  within  the  odd 
{ortho)  rotational  quantum  numbers,  transitions  from  initial  states  other  than  the  ground 
state,  and  transitions  where  the  total  angular  momentum  J  is  not  constrained  to  1/2.  A 
complete  calculation  of  the  S-matrix  elements  over  the  full  distribution  of  states  is 
required  to  obtain  experimentally  useful  information  such  as  cross  sections  and  reaction 
rates. 

Additionally,  one  could  lift  the  restriction  that  P  is  conserved,  thus  ceasing  to  use 
the  centrifugal  sudden  approximation,  which  would  dramatically  increase  the  basis  size. 
A  further  step  would  be  to  allow  for  vibrational  motion  of  the  hydrogen  molecule,  thus 
adding  the  extra  degree  of  freedom  r  back  into  the  problem.  From  this  point,  one  could 
consider  the  full  reactive  problem  B  +  H2  BH  +  H,  with  the  added  difficulty  of  the 
extra  arrangement  channels  permitted  by  the  reaction.  This  requires  a  full  treatment  of 
the  potential  energy  surfaces  for  the  reactive  problem,  as  well  as  significantly  more 
processing  power  to  handle  the  extra  degree  of  freedom. 
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Appendix 

This  appendix  contains  the  full  tabulation  of  probabilities  of  transition  from  the 
ground  state  using  the  three  different  bases.  The  ground  state  has  j  =  0 ,ja  —  1/2,  k  —  0,  and 
0=1/2.  The  probabilities  in  the  ( 0 )  basis,  the  ( 0, 2 )  basis,  and  the  ( 0,  2, 4 )  basis  are 
displayed  in  Figure  21,  Figure  22,  and  Figure  23  respectively. 


Figure  21 .  Probabilities  of  transition  from  the  ground  state  in  the  ( 0 )  basis 
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Figure  22.  Probabilities  of  transition  from  the  ground  state  in  the  ( 0,  2  )  basis 
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Figure  23.  Probabilities  of  transition  from  the  ground  state  in  the  (  0,  2,  4  )  basis 
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